sapply(c(
"rjson",
"data.table",
"dplyr",
"ggplot2",
"stringr",
"purrr",
"foreach",
"patchwork",
"testit",
"lme4",
"lmerTest",
"latex2exp",
"brms"#,
#"tidybayes" # Patchwork/tidybayes ggplot dependency conflict, so source tidybayes directly
),
require, character=TRUE)
## Loading required package: rjson
## Loading required package: data.table
## Loading required package: dplyr
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:data.table':
##
## between, first, last
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
## Loading required package: ggplot2
## Loading required package: stringr
## Loading required package: purrr
##
## Attaching package: 'purrr'
## The following object is masked from 'package:data.table':
##
## transpose
## Loading required package: foreach
##
## Attaching package: 'foreach'
## The following objects are masked from 'package:purrr':
##
## accumulate, when
## Loading required package: patchwork
## Loading required package: testit
## Loading required package: lme4
## Loading required package: Matrix
## Loading required package: lmerTest
##
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
##
## lmer
## The following object is masked from 'package:stats':
##
## step
## Loading required package: latex2exp
## Loading required package: brms
## Loading required package: Rcpp
## Loading 'brms' package (version 2.21.0). Useful instructions
## can be found by typing help('brms'). A more detailed introduction
## to the package is available through vignette('brms_overview').
##
## Attaching package: 'brms'
## The following object is masked from 'package:lme4':
##
## ngrps
## The following object is masked from 'package:stats':
##
## ar
## rjson data.table dplyr ggplot2 stringr purrr foreach
## TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## patchwork testit lme4 lmerTest latex2exp brms
## TRUE TRUE TRUE TRUE TRUE TRUE
sf <- function() invisible(sapply(paste0("./Functions/", list.files("./Functions/", recursive=TRUE)), source)) # Source all fxs
sf()
DefPlotPars()
packageVersion("brms")
## [1] '2.21.0'
Note: These and other code paths need to be adjusted / dir
structure needs to be reconstructed for public data. The public data
also excludes demographics and, due to file size, simulations
(simulations can be rerun using s.R)
Read in learn and test df
learn_df <- read.csv("../data/learn_df.csv")
test_df <- read.csv("../data/test_df.csv")
demogs <- read.csv("../data/demogs_deident.csv")
Demographics
summary(demogs$Age)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 19.00 30.00 37.00 38.42 47.00 65.00
mean(demogs$Age); sd(demogs$Age)
## [1] 38.42029
## [1] 11.42298
table(demogs$Sex)
##
## DATA_EXPIRED Female Male Prefer not to say
## 1 139 133 3
cat("\n Percent female, male, expired, Prefer not to say: \n",
139/nrow(demogs),
133/nrow(demogs),
1/nrow(demogs),
3/nrow(demogs))
##
## Percent female, male, expired, Prefer not to say:
## 0.5036232 0.4818841 0.003623188 0.01086957
table(demogs$Ethnicity.simplified)
##
## Asian Black DATA_EXPIRED Mixed Other White
## 12 26 4 19 12 203
cat("\n % Asian, Black, expired, Mixed, Other, White: \n",
12/nrow(demogs),
26/nrow(demogs),
4/nrow(demogs),
19/nrow(demogs),
12/nrow(demogs),
203/nrow(demogs)
)
##
## % Asian, Black, expired, Mixed, Other, White:
## 0.04347826 0.0942029 0.01449275 0.06884058 0.04347826 0.7355072
Mean delay by set size
learn_df %>% group_by(set_size) %>% summarize(mean(delay, na.rm=TRUE))
## # A tibble: 5 × 2
## set_size `mean(delay, na.rm = TRUE)`
## <int> <dbl>
## 1 1 0
## 2 2 0.832
## 3 3 1.67
## 4 4 2.48
## 5 5 3.33
Load parameter results from best model
## Write into one file and confirm 8/1/24 bug fix was minor
# 11325 and 15031 are the 8/1/24 runs with the minor bug/illogical scaling of epsilon correct
# m35_v1 <- read.csv("../model_res/opt/best/BEST__RunRLWMPRew11325.csv")
# m35_v2 <- read.csv("../model_res/opt/best/BEST__RunRLWMPRew15031.csv")
#ComparePars(m35_v1$nll, m35_v2$nll) # confirm essentially identical
#m35 <- rbind(m35_v1, m35_v2) %>% group_by(ID) %>% slice(which.min(nll))
# And basically identical to older version confirming was a minor bug
# m35_old <- read.csv("../model_res/opt/best/before_8-1-reruns/BEST_m35_RunRLWMPRew.csv")
# ComparePars(m35$nll, m35_old$nll) # confirm essentially identical
#write.csv(m35, "../model_res/opt/best/BEST__m35_RunRLWMPRew-8-1-24-epsilonfixed.csv")
m35 <- read.csv("../model_res/opt/best/BEST__m35_RunRLWMPRew-8-1-24-epsilonfixed.csv") # the final model with the bug/nonsensical implementation of epsilon fixed on 8/1
learn_m <- learn_df %>% group_by(ID) %>% summarize(m=mean(correct))
test_m <- test_df %>% group_by(ID) %>% summarize(m=mean(correct))
assert(all(learn_m$ID==m35$ID))
assert(all(test_m$ID==m35$ID))
And simulations from model
m35_s <-
read.csv("../model_res/sims/SIM_RunRLWMPRew53679.csv")
N simulations
length(unique(m35_s$iter))
## [1] 50
m35_s_learn <- m35_s %>% filter(type=="learning")
m35_s_test <- m35_s %>% filter(type=="test")
pcor_ss <- data.frame(learn_df %>%
group_by(set_size, stim_iter, ID) %>% summarize(m=mean(correct)))
## `summarise()` has grouped output by 'set_size', 'stim_iter'. You can override
## using the `.groups` argument.
pcor_ss_err <- Rmisc::summarySEwithin(pcor_ss,
measurevar = "m",
withinvars = c("set_size", "stim_iter"),
idvar = "ID")
## Automatically converting the following non-factors to factors: set_size, stim_iter
pcor_ss_m <- data.frame(learn_df %>%
group_by(set_size, stim_iter) %>% summarize(m=mean(correct)))
## `summarise()` has grouped output by 'set_size'. You can override using the
## `.groups` argument.
emp_p1 <- ggplot(pcor_ss_err, aes(x=stim_iter, y=m, group=as.factor(set_size), color=as.factor(set_size))) +
geom_line() +
geom_ribbon(aes(ymin=m-se, ymax=m+se), fill='gray57', alpha=.45) +
geom_hline(yintercept=.33, size=1.5, color="gray57") + # chance line
geom_hline(yintercept=c(.5, .6, .7, .8, .9, 1), linetype="dotted") +
geom_vline(xintercept=c(2, 5, 8, 10), linetype="dotted") +
geom_point(aes(fill=as.factor(set_size)), color="black", size=5, pch=21) +
annotate("rect", xmin=6, xmax=10.5, ymin=.3, ymax=1.1, alpha=0.2, fill="gray57") +
ga + ap + lp + xlab("Stimulus iteration") + ylab("Proportion correct") + tol +
ggtitle("Empirical") + tp
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
# Confirm summarySEwithin m and un-normed one are essentially the same
pcor_ss_m <- data.frame(learn_df %>%
group_by(set_size, stim_iter) %>% summarize(m=mean(correct)))
## `summarise()` has grouped output by 'set_size'. You can override using the
## `.groups` argument.
assert(all(round(pcor_ss_err$m, 4)==round(pcor_ss_m$m, 4)))
pcor_ss_sim_m35 <- data.frame(m35_s_learn %>% group_by(set_size, stim_iter) %>%
summarize(m=mean(corrects), n()))
## `summarise()` has grouped output by 'set_size'. You can override using the
## `.groups` argument.
pcor_ss_sim_m35_iter <- data.frame(m35_s_learn %>% filter(iter %in% c(1:30)) %>% group_by(stim_iter, set_size, iter) %>%
summarize(m=mean(corrects), n()))
## `summarise()` has grouped output by 'stim_iter', 'set_size'. You can override
## using the `.groups` argument.
sim_m35_p1 <-
ggplot(pcor_ss_sim_m35, aes(x=as.factor(stim_iter), y=m, group=as.factor(set_size), color=as.factor(set_size))) +
geom_line() +
geom_hline(yintercept=.33, size=1.5, color="gray57") + # chance line
geom_hline(yintercept=c(.5, .6, .7, .8, .9, 1), linetype="dotted") +
geom_vline(xintercept=c(2, 5, 8, 10), linetype="dotted") +
geom_jitter(data=pcor_ss_sim_m35_iter, aes(fill=as.factor(set_size)), color="black", height=0, width=.2, alpha=1, size=2, pch=21) +
geom_point(aes(fill=as.factor(set_size)), color="black", size=6, pch=21, alpha=.7) +
annotate("rect", xmin=6, xmax=10.5, ymin=.3, ymax=1.1, alpha=0.2, fill="gray57") +
ga + ap + lp + xlab("Stimulus iteration") + ylab("") +
tp + ggtitle("Simulated")
emp_sim_perf <- emp_p1 + sim_m35_p1
emp_sim_perf
#ggsave("../paper/figs/pieces/fig2_emp_sim_perf.png", emp_sim_perf, height = 3.5, width=11, dpi=300)
Before/after break regressor
end_p1_begin_p2_df <- learn_df %>% filter(stim_iter %in% c(5, 6))
#end_p1_begin_p2_df$before_after_break <- end_p1_begin_p2_df$stim_iter
end_p1_begin_p2_df[end_p1_begin_p2_df$stim_iter==5, "before_after_break"] <- -1
end_p1_begin_p2_df[end_p1_begin_p2_df$stim_iter==6, "before_after_break"] <- 1
Main effects of set size and stim iter and interaction reflecting parametric effect over time in both phase 1 and phase 2
summary(pcor_phase_1_v2 <-
glmer(correct ~ scale(set_size)*scale(stim_iter) +
(scale(set_size)*scale(stim_iter)|ID),
data=learn_df %>% filter(phase==1),
family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: correct ~ scale(set_size) * scale(stim_iter) + (scale(set_size) *
## scale(stim_iter) | ID)
## Data: learn_df %>% filter(phase == 1)
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 44084.9 44205.7 -22028.5 44056.9 41236
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -11.8169 -0.7807 0.3279 0.6811 1.7107
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 0.19502 0.4416
## scale(set_size) 0.04951 0.2225 -0.08
## scale(stim_iter) 0.10059 0.3172 0.90 -0.23
## scale(set_size):scale(stim_iter) 0.02250 0.1500 -0.16 0.69 -0.53
## Number of obs: 41250, groups: ID, 275
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.99731 0.03038 32.83 <2e-16 ***
## scale(set_size) -0.54013 0.02075 -26.03 <2e-16 ***
## scale(stim_iter) 1.07075 0.02424 44.18 <2e-16 ***
## scale(set_size):scale(stim_iter) -0.37241 0.01824 -20.42 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) scl(st_s) scl(stm_)
## scal(st_sz) -0.201
## scl(stm_tr) 0.774 -0.297
## scl(s_):(_) -0.228 0.587 -0.417
summary(pcor_phase_2_v2 <-
glmer(correct ~ scale(set_size)*scale(stim_iter) +
(scale(set_size)*scale(stim_iter)|ID),
data=learn_df %>% filter(phase==2),
family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: correct ~ scale(set_size) * scale(stim_iter) + (scale(set_size) *
## scale(stim_iter) | ID)
## Data: learn_df %>% filter(phase == 2)
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 32506.6 32627.4 -16239.3 32478.6 41236
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -9.0386 0.1654 0.3001 0.4631 1.3998
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 0.61266 0.7827
## scale(set_size) 0.11233 0.3352 0.57
## scale(stim_iter) 0.02303 0.1517 0.56 0.36
## scale(set_size):scale(stim_iter) 0.01768 0.1330 0.26 0.24 -0.56
## Number of obs: 41250, groups: ID, 275
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.02297 0.05113 39.568 <2e-16 ***
## scale(set_size) -0.28678 0.02877 -9.969 <2e-16 ***
## scale(stim_iter) 0.69938 0.02066 33.849 <2e-16 ***
## scale(set_size):scale(stim_iter) -0.25844 0.02082 -12.410 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) scl(st_s) scl(stm_)
## scal(st_sz) 0.298
## scl(stm_tr) 0.410 -0.049
## scl(s_):(_) 0.003 0.436 -0.291
# Singular
# summary(pcor_before_after <-
# glmer(correct ~ scale(set_size)*scale(before_after_break) +
# (scale(set_size)*scale(before_after_break)|ID),
# data=end_p1_begin_p2_df,
# family="binomial", control = glmerControl(optimizer = "bobyqa")))
summary(pcor_before_after <-
glmer(correct ~ scale(set_size)*scale(before_after_break) +
(scale(set_size) + scale(before_after_break)|ID),
data=end_p1_begin_p2_df,
family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## correct ~ scale(set_size) * scale(before_after_break) + (scale(set_size) +
## scale(before_after_break) | ID)
## Data: end_p1_begin_p2_df
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 16446.4 16523.6 -8213.2 16426.4 16490
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.7814 0.1644 0.3936 0.5720 1.5714
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 0.421477 0.64921
## scale(set_size) 0.068597 0.26191 0.31
## scale(before_after_break) 0.002658 0.05156 0.34 -0.44
## Number of obs: 16500, groups: ID, 275
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.41359 0.04558 31.02 <2e-16
## scale(set_size) -0.29171 0.02908 -10.03 <2e-16
## scale(before_after_break) -0.52282 0.02316 -22.58 <2e-16
## scale(set_size):scale(before_after_break) 0.37378 0.02349 15.91 <2e-16
##
## (Intercept) ***
## scale(set_size) ***
## scale(before_after_break) ***
## scale(set_size):scale(before_after_break) ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) scl(_) sc(__)
## scal(st_sz) 0.025
## scl(bfr_f_) -0.146 0.229
## scl(_):(__) 0.156 -0.362 -0.277
car::vif(pcor_before_after)
## scale(set_size)
## 1.174968
## scale(before_after_break)
## 1.105577
## scale(set_size):scale(before_after_break)
## 1.205734
car::vif(pcor_phase_2_v2)
## scale(set_size) scale(stim_iter)
## 1.245388 1.101492
## scale(set_size):scale(stim_iter)
## 1.357267
car::vif(pcor_phase_1_v2)
## scale(set_size) scale(stim_iter)
## 1.533658 1.216770
## scale(set_size):scale(stim_iter)
## 1.692973
Logistic model predictions are a bit mis-specified but decent
sjPlot::plot_model(pcor_phase_1_v2,
type = "pred", terms = c("stim_iter", "set_size")) +
ga + ap + lp + tp + xlab("Stimulus iteration") +
ylab("Proportion correct") + ggtitle("Phase 1 regression model predictions")
sjPlot::tab_model(pcor_phase_1_v2)
| correct | |||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 2.71 | 2.55 – 2.88 | <0.001 |
| set size | 0.58 | 0.56 – 0.61 | <0.001 |
| stim iter | 2.92 | 2.78 – 3.06 | <0.001 |
| set size × stim iter | 0.69 | 0.66 – 0.71 | <0.001 |
| Random Effects | |||
| σ2 | 3.29 | ||
| τ00 ID | 0.20 | ||
| τ11 ID.scale(set_size) | 0.05 | ||
| τ11 ID.scale(stim_iter) | 0.10 | ||
| τ11 ID.scale(set_size):scale(stim_iter) | 0.02 | ||
| ρ01 | -0.08 | ||
| 0.90 | |||
| -0.16 | |||
| ICC | 0.10 | ||
| N ID | 275 | ||
| Observations | 41250 | ||
| Marginal R2 / Conditional R2 | 0.301 / 0.371 | ||
sjPlot::plot_model(pcor_phase_2_v2,
type = "pred", terms = c("stim_iter", "set_size")) +
ga + ap + lp + tp + xlab("Stimulus iteration") +
ylab("Proportion correct") + ggtitle("Phase 2 regression model predictions")
sjPlot::tab_model(pcor_phase_2_v2)
| correct | |||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 7.56 | 6.84 – 8.36 | <0.001 |
| set size | 0.75 | 0.71 – 0.79 | <0.001 |
| stim iter | 2.01 | 1.93 – 2.10 | <0.001 |
| set size × stim iter | 0.77 | 0.74 – 0.80 | <0.001 |
| Random Effects | |||
| σ2 | 3.29 | ||
| τ00 ID | 0.61 | ||
| τ11 ID.scale(set_size) | 0.11 | ||
| τ11 ID.scale(stim_iter) | 0.02 | ||
| τ11 ID.scale(set_size):scale(stim_iter) | 0.02 | ||
| ρ01 | 0.57 | ||
| 0.56 | |||
| 0.26 | |||
| ICC | 0.19 | ||
| N ID | 275 | ||
| Observations | 41250 | ||
| Marginal R2 / Conditional R2 | 0.136 / 0.299 | ||
pcor_before_after_p <- sjPlot::plot_model(pcor_before_after,
type = "pred", terms = c("before_after_break", "set_size")) +
ga + ap + lp + tp + xlab("Stimulus iteration") +
ylab("Proportion correct") + ggtitle("Before to after break regression model predictions")
pcor_before_after_p
sjPlot::tab_model(pcor_before_after)
| correct | |||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 4.11 | 3.76 – 4.49 | <0.001 |
| set size | 0.75 | 0.71 – 0.79 | <0.001 |
| before after break | 0.59 | 0.57 – 0.62 | <0.001 |
|
set size × before after break |
1.45 | 1.39 – 1.52 | <0.001 |
| Random Effects | |||
| σ2 | 3.29 | ||
| τ00 ID | 0.42 | ||
| τ11 ID.scale(set_size) | 0.07 | ||
| τ11 ID.scale(before_after_break) | 0.00 | ||
| ρ01 | 0.31 | ||
| 0.34 | |||
| ICC | 0.13 | ||
| N ID | 275 | ||
| Observations | 16500 | ||
| Marginal R2 / Conditional R2 | 0.116 / 0.231 | ||
#ggsave("../paper/figs/supp-figs/pcor_before_after_p.png", pcor_before_after_p, height = 5, width=11, dpi=300)
Trial 1 so there’s def no S-A-O load on WM
rt_trial_1 <- data.frame(learn_df %>% filter(stim_iter==1) %>%
group_by(set_size, ID) %>% summarize(m=mean(rt)))
## `summarise()` has grouped output by 'set_size'. You can override using the
## `.groups` argument.
rt_trial_1_m <- data.frame(learn_df %>% filter(stim_iter==1) %>%
group_by(set_size) %>% summarize(m=mean(rt)))
rt_tr1_err <- Rmisc::summarySEwithin(rt_trial_1,
measurevar = "m",
withinvars = c("set_size"),
idvar = "ID")
## Automatically converting the following non-factors to factors: set_size
assert(all(round(rt_tr1_err$m, 4)==round(rt_trial_1_m$m, 4)))
emp_rt_tr_1 <- ggplot(rt_tr1_err ,
aes(x=set_size, y=m, group=as.factor(set_size), fill=as.factor(set_size))) +
geom_errorbar(aes(x=set_size, ymin=m-se, ymax=m+se), size=1.5, width=.25) +
geom_point(size=6, pch=21, color="black", alpha=.7) +
ga + ap + lp + xlab("Set size") + ylab("") +
ggtitle("At trial 1") + tp + tol#+ ylim(600, 900)
rt_ss_si6 <- data.frame(learn_df %>% filter(stim_iter==6) %>%
group_by(set_size, ID) %>% summarize(m=mean(rt)))
## `summarise()` has grouped output by 'set_size'. You can override using the
## `.groups` argument.
rt_ss_si6_m <- data.frame(learn_df %>% filter(stim_iter==6) %>%
group_by(set_size) %>% summarize(m=mean(rt)))
rt_ss_si6_err <- Rmisc::summarySEwithin(rt_ss_si6,
measurevar = "m",
withinvars = c("set_size"),
idvar = "ID")
## Automatically converting the following non-factors to factors: set_size
assert(all(round(rt_ss_si6_err$m, 4)==round(rt_ss_si6_m$m, 4)))
emp_rt_si_6 <- ggplot(rt_ss_si6_err ,
aes(x=set_size, y=m, group=as.factor(set_size), fill=as.factor(set_size))) +
geom_errorbar(aes(x=set_size, ymin=m-se, ymax=m+se), size=1.5, width=.25) +
geom_point(size=6, pch=21, color="black", alpha=.7) +
ga + ap + lp + xlab("Set size") + ylab("") +
ggtitle("At stimulus iteration 6") + tp + tol
rt_ss_si6_df <- data.frame(learn_df %>% filter(stim_iter==6))
rt_ss_si6_df$set_size_indicator <- if_else(rt_ss_si6_df$set_size==1, 1, -1)
# Reduced because model without RE did not converge
summary(rt_ss_si6_df_model <-
lmer(rt ~ set_size_indicator + (1|ID), data=rt_ss_si6_df))
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: rt ~ set_size_indicator + (1 | ID)
## Data: rt_ss_si6_df
##
## REML criterion at convergence: 111694.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.5966 -0.6484 -0.0592 0.6007 3.7289
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 7416 86.11
## Residual 41809 204.47
## Number of obs: 8250, groups: ID, 275
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 737.512 6.880 593.431 107.20 <2e-16 ***
## set_size_indicator 49.475 4.512 7974.000 10.96 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## st_sz_ndctr 0.568
rt_ss <- data.frame(learn_df %>% group_by(set_size, stim_iter, ID) %>% summarize(m=mean(rt)))
## `summarise()` has grouped output by 'set_size', 'stim_iter'. You can override
## using the `.groups` argument.
rt_ss_m <- data.frame(learn_df %>%
group_by(set_size, stim_iter) %>% summarize(m=mean(rt)))
## `summarise()` has grouped output by 'set_size'. You can override using the
## `.groups` argument.
rt_ss_err <- Rmisc::summarySEwithin(rt_ss,
measurevar = "m",
withinvars = c("set_size", "stim_iter"),
idvar = "ID")
## Automatically converting the following non-factors to factors: set_size, stim_iter
assert(all(round(rt_ss_err$m, 4)==round(rt_ss_m$m, 4)))
emp_rt <-
ggplot(rt_ss_err, aes(x=stim_iter, y=m, group=as.factor(set_size), color=as.factor(set_size))) +
geom_line() +
geom_ribbon(aes(ymin=m-se, ymax=m+se), fill='gray57', alpha=.45) +
geom_vline(xintercept=c(3, 7, 9), linetype="dotted") +
geom_point(aes(fill=as.factor(set_size)), color="black", size=5, pch=21) +
annotate("rect", xmin=6, xmax=10.5, ymin=300, ymax=850, alpha=0.2, fill="gray57") +
ga + ap + lp + xlab("Stimulus iteration") + ylab("Reaction time") + tol #+
#ggtitle("") + tp #+ ylim(300, 900)
rt_ss_d0_r1 <-
data.frame(learn_df %>% filter(delay==0 & past_reward == 1) %>%
group_by(set_size, stim_iter, ID) %>% summarize(m=mean(rt)))
## `summarise()` has grouped output by 'set_size', 'stim_iter'. You can override
## using the `.groups` argument.
rt_ss_d0_r1_m <- data.frame(learn_df %>% filter(delay==0 & past_reward == 1) %>%
group_by(set_size, stim_iter) %>% summarize(m=mean(rt)))
## `summarise()` has grouped output by 'set_size'. You can override using the
## `.groups` argument.
rt_ss_err_d0_r1 <- Rmisc::summarySEwithin(rt_ss_d0_r1,
measurevar = "m",
withinvars = c("set_size", "stim_iter"),
idvar = "ID")
## Automatically converting the following non-factors to factors: set_size, stim_iter
assert(all(round(rt_ss_err_d0_r1$m, 4)==round(rt_ss_err_d0_r1$m, 4)))
emp_rt_d0_r1 <- ggplot(rt_ss_err_d0_r1, aes(x=stim_iter, y=m, group=as.factor(set_size), color=as.factor(set_size))) +
geom_line() +
geom_ribbon(aes(ymin=m-se, ymax=m+se), fill='gray57', alpha=.45) +
geom_vline(xintercept=c(2, 5, 7), linetype="dotted") +
geom_point(aes(fill=as.factor(set_size)), color="black", size=5, pch=21) +
annotate("rect", xmin=5, xmax=8.5, ymin=300, ymax=630, alpha=0.2, fill="gray57") +
ga + ap + lp + xlab("Stimulus iteration") + ylab("") +
ggtitle("When prior reward and no delay") + tp + tol #+ ylim(300, 900)
rt_nodrew <- learn_df %>% filter(delay==0 & past_reward == 1)
summary(rt_nodrew_model <-
lmer(rt ~ scale(set_size) + (scale(set_size)|ID), data=rt_nodrew))
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: rt ~ scale(set_size) + (scale(set_size) | ID)
## Data: rt_nodrew
##
## REML criterion at convergence: 239851.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.1271 -0.6197 -0.1287 0.4544 6.7795
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 4851.9 69.66
## scale(set_size) 480.4 21.92 -0.08
## Residual 20929.7 144.67
## Number of obs: 18678, groups: ID, 275
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 480.237 4.334 274.098 110.81 <2e-16 ***
## scale(set_size) 67.556 1.699 271.451 39.76 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## scal(st_sz) -0.059
all_rt <-
emp_rt + emp_rt_d0_r1 + emp_rt_tr_1 + emp_rt_si_6 + plot_annotation(title="Empirical reaction time", theme = theme(plot.title = element_text(size = 25, hjust=.5)))#, size=10)
all_rt
#ggsave("../paper/figs/pieces/fig2_emp_rt.png", all_rt, height = 6, width=11, dpi=300)
Confirms strong set size effect even at first stim iter — suggesting different proactive strategies
summary(rt_si_1 <-
lmer(rt ~ scale(set_size) + (scale(set_size)|ID), data=learn_df %>% filter(stim_iter==1)))
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: rt ~ scale(set_size) + (scale(set_size) | ID)
## Data: learn_df %>% filter(stim_iter == 1)
##
## REML criterion at convergence: 112275.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.1250 -0.6697 -0.0766 0.5833 4.1575
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 9725 98.62
## scale(set_size) 1063 32.61 0.31
## Residual 43818 209.33
## Number of obs: 8250, groups: ID, 275
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 716.300 6.378 274.026 112.3 <2e-16 ***
## scale(set_size) 37.568 3.030 273.996 12.4 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## scal(st_sz) 0.188
sjPlot::tab_model(rt_si_1)
| rt | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 716.30 | 703.80 – 728.80 | <0.001 |
| set size | 37.57 | 31.63 – 43.51 | <0.001 |
| Random Effects | |||
| σ2 | 43818.07 | ||
| τ00 ID | 9725.25 | ||
| τ11 ID.scale(set_size) | 1063.36 | ||
| ρ01 ID | 0.31 | ||
| ICC | 0.20 | ||
| N ID | 275 | ||
| Observations | 8250 | ||
| Marginal R2 / Conditional R2 | 0.025 / 0.218 | ||
Regression models in phase 1 and phase 2 show main effects of set size and stim iter and an interaction in both phases
(REs complexity reduce due to convergence failure with full)
summary(rt_p1 <-
lmer(rt ~ scale(set_size)*scale(stim_iter) +
(scale(set_size) + scale(stim_iter)|ID),
data=learn_df %>% filter(phase==1)))
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: rt ~ scale(set_size) * scale(stim_iter) + (scale(set_size) +
## scale(stim_iter) | ID)
## Data: learn_df %>% filter(phase == 1)
##
## REML criterion at convergence: 558890.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.0355 -0.6793 -0.1269 0.5603 4.7845
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 6923.2 83.21
## scale(set_size) 951.4 30.84 0.33
## scale(stim_iter) 428.5 20.70 -0.20 -0.18
## Residual 43257.6 207.98
## Number of obs: 41250, groups: ID, 275
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 648.359 5.121 274.006 126.61
## scale(set_size) 88.066 2.123 274.000 41.48
## scale(stim_iter) -38.968 1.615 274.014 -24.14
## scale(set_size):scale(stim_iter) 18.135 1.024 40424.000 17.71
## Pr(>|t|)
## (Intercept) <2e-16 ***
## scale(set_size) <2e-16 ***
## scale(stim_iter) <2e-16 ***
## scale(set_size):scale(stim_iter) <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) scl(st_s) scl(stm_)
## scal(st_sz) 0.286
## scl(stm_tr) -0.154 -0.124
## scl(s_):(_) 0.000 0.000 0.000
sjPlot::tab_model(rt_p1)
| rt | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 648.36 | 638.32 – 658.40 | <0.001 |
| set size | 88.07 | 83.90 – 92.23 | <0.001 |
| stim iter | -38.97 | -42.13 – -35.80 | <0.001 |
| set size × stim iter | 18.13 | 16.13 – 20.14 | <0.001 |
| Random Effects | |||
| σ2 | 43257.63 | ||
| τ00 ID | 6923.16 | ||
| τ11 ID.scale(set_size) | 951.39 | ||
| τ11 ID.scale(stim_iter) | 428.50 | ||
| ρ01 | 0.33 | ||
| -0.20 | |||
| ICC | 0.16 | ||
| N ID | 275 | ||
| Observations | 41250 | ||
| Marginal R2 / Conditional R2 | 0.157 / 0.293 | ||
sjPlot::plot_model(rt_p1, type = "pred", terms = c("stim_iter", "set_size")) +
ga + ap + lp + tp + xlab("Stimulus iteration") + ylab("Reaction time") + ggtitle("Phase 1 regression model predictions")
summary(rt_p2 <-lmer(rt ~ scale(set_size)*scale(stim_iter) + (scale(set_size) + scale(stim_iter)|ID), data=learn_df %>% filter(phase==2)))
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: rt ~ scale(set_size) * scale(stim_iter) + (scale(set_size) +
## scale(stim_iter) | ID)
## Data: learn_df %>% filter(phase == 2)
##
## REML criterion at convergence: 545785.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.5886 -0.6438 -0.1083 0.5057 7.5230
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 6462.4 80.39
## scale(set_size) 626.8 25.04 0.28
## scale(stim_iter) 143.2 11.97 -0.12 -0.19
## Residual 31520.5 177.54
## Number of obs: 41250, groups: ID, 275
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 580.8416 4.9258 273.9913 117.92
## scale(set_size) 57.9437 1.7445 274.0246 33.22
## scale(stim_iter) -40.3142 1.1336 273.9947 -35.56
## scale(set_size):scale(stim_iter) 19.6083 0.8742 40423.9930 22.43
## Pr(>|t|)
## (Intercept) <2e-16 ***
## scale(set_size) <2e-16 ***
## scale(stim_iter) <2e-16 ***
## scale(set_size):scale(stim_iter) <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) scl(st_s) scl(stm_)
## scal(st_sz) 0.234
## scl(stm_tr) -0.075 -0.104
## scl(s_):(_) 0.000 0.000 0.000
car::vif(rt_p1)
## scale(set_size) scale(stim_iter)
## 1.015515 1.015515
## scale(set_size):scale(stim_iter)
## 1.000000
car::vif(rt_p2)
## scale(set_size) scale(stim_iter)
## 1.010985 1.010985
## scale(set_size):scale(stim_iter)
## 1.000000
sjPlot::tab_model(rt_p2)
| rt | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 580.84 | 571.19 – 590.50 | <0.001 |
| set size | 57.94 | 54.52 – 61.36 | <0.001 |
| stim iter | -40.31 | -42.54 – -38.09 | <0.001 |
| set size × stim iter | 19.61 | 17.89 – 21.32 | <0.001 |
| Random Effects | |||
| σ2 | 31520.54 | ||
| τ00 ID | 6462.37 | ||
| τ11 ID.scale(set_size) | 626.77 | ||
| τ11 ID.scale(stim_iter) | 143.23 | ||
| ρ01 | 0.28 | ||
| -0.12 | |||
| ICC | 0.19 | ||
| N ID | 275 | ||
| Observations | 41250 | ||
| Marginal R2 / Conditional R2 | 0.122 / 0.286 | ||
sjPlot::plot_model(rt_p2, type = "pred", terms = c("stim_iter", "set_size")) +
ga + ap + lp + tp + xlab("Stimulus iteration") + ylab("Reaction time") + ggtitle("Phase 2 regression model predictions")
learn_ms <- data.frame(learn_df %>%
group_by(ID, set_size) %>% summarize(learn_m=mean(correct)))
## `summarise()` has grouped output by 'ID'. You can override using the `.groups`
## argument.
test_ms <- data.frame(test_df %>%
group_by(ID, set_size) %>% summarize(test_m=mean(correct)))
## `summarise()` has grouped output by 'ID'. You can override using the `.groups`
## argument.
assert(all(test_ms$ID==learn_ms$ID))
test_ms$learn_m <- learn_ms$learn_m
emp_corrs <- ggplot(test_ms,
aes(x=learn_m, y=test_m, color=as.factor(set_size))) +
geom_jitter(alpha=.4, size=2, width=.02) +
ggpubr::stat_cor(method="spearman", size=5, label.y=1.12) + ga + ft + ap +
xlab("Learning") + ylab("Test") + facet_wrap(~ set_size, nrow=1) + tol +
xlim(0, 1) + ylim(0, 1.17) +
theme(axis.text = element_blank(), axis.ticks = element_blank()) +
ggtitle("Proportion correct correlation \n Empirical") + tp
Take one simulation iter to match number used empirically and not “smooth” by more sims
learn_ms_s <- data.frame(m35_s_learn %>% filter(iter==1) %>% group_by(ID, set_size) %>% summarize(learn_m=mean(corrects)))
## `summarise()` has grouped output by 'ID'. You can override using the `.groups`
## argument.
test_ms_s <- data.frame(m35_s_test %>% filter(iter==1) %>% group_by(ID, set_size) %>% summarize(test_m=mean(corrects)))
## `summarise()` has grouped output by 'ID'. You can override using the `.groups`
## argument.
assert(all(test_ms_s$ID==learn_ms_s$ID))
test_ms_s$learn_m <- learn_ms_s$learn_m
sim_corrs <- ggplot(test_ms_s, aes(x=learn_m, y=test_m, color=as.factor(set_size))) + geom_jitter(alpha=.4, size=2, width=.0) +
ggpubr::stat_cor(method="spearman", size=5, label.y=1.12) + ga + ft + ap +
xlab("Learning") + ylab("Test") + facet_wrap(~ set_size, nrow=1) + tol +
xlim(0, 1) + ylim(0, 1.17) +
theme(axis.text = element_blank(), axis.ticks = element_blank()) +
ggtitle("Simulated") + tp
A few missing due to jitter and divide by 0
all_corrs <- emp_corrs / sim_corrs
all_corrs
## Warning: Removed 17 rows containing missing values (`geom_point()`).
#ggsave("../paper/figs/pieces/fig2_perf_corrs.png", all_corrs, height = 6, width=11, dpi=300)
Empirical
p1_test_trials <- test_df %>% filter(phase==1)
p2_test_trials <- test_df %>% filter(phase==2)
pcor_p1_test <-
data.frame(p1_test_trials %>% group_by(set_size, ID) %>% summarize(m=mean(correct)))
## `summarise()` has grouped output by 'set_size'. You can override using the
## `.groups` argument.
pcor_p2_test <-
data.frame(p2_test_trials %>% group_by(set_size, ID) %>% summarize(m=mean(correct)))
## `summarise()` has grouped output by 'set_size'. You can override using the
## `.groups` argument.
pcor_si6 <- learn_df %>% filter(stim_iter==6) %>%
group_by(set_size, ID) %>% summarize(m=mean(correct))
## `summarise()` has grouped output by 'set_size'. You can override using the
## `.groups` argument.
pcor_p1_test_err <- Rmisc::summarySEwithin(pcor_p1_test,
measurevar = "m",
withinvars = c("set_size"),
idvar = "ID")
## Automatically converting the following non-factors to factors: set_size
pcor_si6_err <- Rmisc::summarySEwithin(pcor_si6,
measurevar = "m",
withinvars = c("set_size"),
idvar = "ID")
## Automatically converting the following non-factors to factors: set_size
pcor_p2_test_err <- Rmisc::summarySEwithin(pcor_p2_test,
measurevar = "m",
withinvars = c("set_size"),
idvar = "ID")
## Automatically converting the following non-factors to factors: set_size
# Continuing to sanity check summarSEwithin
pcor_p1_test_m <-
data.frame(p1_test_trials %>% group_by(set_size) %>% summarize(m=mean(correct)))
pcor_p2_test_m <-
data.frame(p2_test_trials %>% group_by(set_size) %>% summarize(m=mean(correct)))
pcor_si6_m <- learn_df %>% filter(stim_iter==6) %>%
group_by(set_size) %>% summarize(m=mean(correct))
assert(all(round(pcor_p1_test_err$m, 4)==round(pcor_p1_test_m$m, 4)))
assert(all(round(pcor_p2_test_err$m, 4)==round(pcor_p2_test_m$m, 4)))
assert(all(round(pcor_si6_err$m, 4)==round(pcor_si6_m$m, 4)))
emp_p1_test <- ggplot(pcor_p1_test_err, aes(x=set_size, y=m, fill=set_size)) +
geom_hline(yintercept = seq(.1, 1, .1), alpha=.3) +
geom_bar(stat="identity", color="black") +
geom_errorbar(aes(ymin=m-se, ymax=m+se), width=.2) +
ga + ap + xlab("Set size") + ylab("Proportion correct") + tol + ylim(0, 1) +
tp + ggtitle("Test phase 1")
emp_si6 <- ggplot(pcor_si6_err, aes(x=set_size, y=m, fill=set_size)) +
geom_hline(yintercept = seq(.1, 1, .1), alpha=.3) +
geom_bar(stat="identity", color="black") +
geom_errorbar(aes(ymin=m-se, ymax=m+se), width=.2) +
ga + ap + xlab("Set size") + ylab("") + tol + ylim(0, 1) +
tp + ggtitle("Stimulus iteration 6")
emp_p2_test <- ggplot(pcor_p2_test_err, aes(x=set_size, y=m, fill=set_size)) +
geom_hline(yintercept = seq(.1, 1, .1), alpha=.3) +
geom_bar(stat="identity", color="black") +
geom_errorbar(aes(ymin=m-se, ymax=m+se), width=.2) +
ga + ap + xlab("Set size") + ylab("") + tol + ylim(0, 1) + tp +
ggtitle("Test phase 2")
U-shaped pattern evident at test phase 1, which replicates at the beginning of the next learning phase (stimulus iteration 6) before there has been a chance for further learning (thus making it effectively another test trial)
The effect then appears to reduce substantially by test phase 2 — with stim iter 1 still decreased, but set size 2 if anything higher than the others
emps_u_plot <-
emp_p1_test + emp_si6 + emp_p2_test + plot_annotation(title="Empirical", theme = theme(plot.title = element_text(size = 25, hjust=.5)))#,
emps_u_plot
Simulation plots
pcor_p1_test_sim_m35 <-
data.frame(m35_s_test %>% filter(phase==1) %>% group_by(set_size) %>% summarize(m=mean(corrects)))
pcor_p1_test_sim_m35_iters <-
data.frame(m35_s_test %>% filter(phase==1) %>% group_by(set_size, iter) %>% summarize(m=mean(corrects)))
## `summarise()` has grouped output by 'set_size'. You can override using the
## `.groups` argument.
pcor_si6_sim_m35 <-
data.frame(m35_s_learn %>% filter(stim_iter==6) %>% group_by(set_size) %>% summarize(m=mean(corrects)))
pcor_si6_sim_m35_iters <-
data.frame(m35_s_learn %>% filter(stim_iter==6) %>% group_by(set_size, iter) %>% summarize(m=mean(corrects)))
## `summarise()` has grouped output by 'set_size'. You can override using the
## `.groups` argument.
pcor_p2_test_sim_m35 <-
data.frame(m35_s_test %>% filter(phase==2) %>% group_by(set_size) %>% summarize(m=mean(corrects)))
pcor_p2_test_sim_m35_iters <-
data.frame(m35_s_test %>% filter(phase==2) %>% group_by(set_size, iter) %>% summarize(m=mean(corrects)))
## `summarise()` has grouped output by 'set_size'. You can override using the
## `.groups` argument.
sim_p1_test_m35 <-
ggplot(pcor_p1_test_sim_m35, aes(x=as.factor(set_size), y=m, fill=as.factor(set_size))) +
geom_hline(yintercept = seq(.1, 1, .1), alpha=.3) +
geom_bar(stat="identity", color="black") +
geom_jitter(data=pcor_p1_test_sim_m35_iters, size=2, alpha=1, width=.08, height=0, pch=21,
aes(x=as.factor(set_size), y=m, fill=as.factor(set_size))) +
ga + ap + xlab("Set size") + ylab("Proportion correct") + tol + ylim(0, 1) +
tp + ggtitle("Test phase 1")
sim_si6_m35 <- ggplot(pcor_si6_sim_m35, aes(x=as.factor(set_size), y=m, fill=as.factor(set_size))) +
geom_hline(yintercept = seq(.1, 1, .1), alpha=.3) +
geom_bar(stat="identity", color="black") +
geom_jitter(data=pcor_si6_sim_m35_iters, size=2, alpha=1, width=.08, height=0, pch=21,
aes(x=as.factor(set_size), y=m, fill=as.factor(set_size))) +
ga + ap + xlab("Set size") + ylab("") + tol + ylim(0, 1) +
tp + ggtitle("Stimulus iteration 6")
sim_p2_test_m35 <- ggplot(pcor_p2_test_sim_m35, aes(x=as.factor(set_size), y=m, fill=as.factor(set_size))) +
geom_hline(yintercept = seq(.1, 1, .1), alpha=.3) +
geom_bar(stat="identity", color="black") +
geom_jitter(data=pcor_p2_test_sim_m35_iters, size=2, alpha=1, width=.08, height=0, pch=21,
aes(x=as.factor(set_size), y=m, fill=as.factor(set_size))) +
ga + ap + xlab("Set size") + ylab("") + tol + ylim(0, 1) + tp +
tp + ggtitle("Test phase 2")
sims_u_plot <-
sim_p1_test_m35 + sim_si6_m35 + sim_p2_test_m35 + plot_annotation(title="Simulated", theme = theme(plot.title = element_text(size = 25, hjust=.5)))#,
sims_u_plot
# ggsave("../paper/figs/pieces/fig3_empU.png", emps_u_plot, height = 4, width=11, dpi=300)
#ggsave("../paper/figs/pieces/fig3_simU.png", sims_u_plot, height = 4, width=11, dpi=300)
summary(pcor_p1_test_quad <-
glmer(correct ~ poly(set_size, 2) +
(poly(set_size, 2)|ID),
data=test_df %>% filter(phase==1),
family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: correct ~ poly(set_size, 2) + (poly(set_size, 2) | ID)
## Data: test_df %>% filter(phase == 1)
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 18853.1 18922.5 -9417.5 18835.1 16491
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.9057 -0.8769 0.4293 0.6433 1.9108
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 0.6646 0.8152
## poly(set_size, 2)1 1999.3502 44.7141 0.15
## poly(set_size, 2)2 2394.8985 48.9377 -0.13 -0.12
## Number of obs: 16500, groups: ID, 275
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.95321 0.05283 18.043 < 2e-16 ***
## poly(set_size, 2)1 13.42438 3.24378 4.138 3.50e-05 ***
## poly(set_size, 2)2 -26.23676 3.43764 -7.632 2.31e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) p(_,2)1
## ply(st_,2)1 0.096
## ply(st_,2)2 -0.099 -0.065
Predictions look reasonable and do fit a clear inverted U
sjPlot::tab_model(pcor_p1_test_quad)
| correct | |||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 2.59 | 2.34 – 2.88 | <0.001 |
| set size [1st degree] | 676291.41 | 1172.25 – 390163629.22 | <0.001 |
| set size [2nd degree] | 0.00 | 0.00 – 0.00 | <0.001 |
| Random Effects | |||
| σ2 | 3.29 | ||
| τ00 ID | 0.66 | ||
| τ11 ID.poly(set_size, 2)1 | 1999.35 | ||
| τ11 ID.poly(set_size, 2)2 | 2394.90 | ||
| ρ01 | 0.15 | ||
| -0.13 | |||
| ICC | 0.22 | ||
| N ID | 275 | ||
| Observations | 16500 | ||
| Marginal R2 / Conditional R2 | 0.012 / 0.230 | ||
pcor_p1_test_quad_p <- sjPlot::plot_model(pcor_p1_test_quad, type = "pred",
terms = c("set_size [all]")) +
ga + ap + lp + tp + xlab("Set size") +
ylab("Proportion correct") + ggtitle("Quadratic model regression predictions \n Phase 1 test")
pcor_p1_test_quad_p
Stimulus iteration 6 replicates effect
summary(pcor_si6_quad <-
glmer(correct ~ poly(set_size, 2) +
(poly(set_size, 2)|ID),
data=learn_df %>% filter(stim_iter==6),
family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: correct ~ poly(set_size, 2) + (poly(set_size, 2) | ID)
## Data: learn_df %>% filter(stim_iter == 6)
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 9777.1 9840.3 -4879.6 9759.1 8241
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9925 -1.0132 0.4758 0.6605 1.4247
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 0.4794 0.6924
## poly(set_size, 2)1 414.1219 20.3500 0.28
## poly(set_size, 2)2 366.8434 19.1532 -0.48 0.62
## Number of obs: 8250, groups: ID, 275
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.91428 0.04948 18.478 < 2e-16 ***
## poly(set_size, 2)1 6.44043 2.43970 2.640 0.00829 **
## poly(set_size, 2)2 -22.80614 2.48653 -9.172 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) p(_,2)1
## ply(st_,2)1 0.122
## ply(st_,2)2 -0.237 0.204
pcor_si6_quad_p <- sjPlot::plot_model(pcor_si6_quad,
type = "pred", terms = c("set_size")) +
ga + ap + lp + tp + xlab("Set size") +
ylab("Proportion correct") +
ggtitle("Quadratic model regression predictions \n Stim iter 6")
## Model contains splines or polynomial terms. Consider using
## `terms="set_size [all]"` to get smooth plots. See also package-vignette
## 'Marginal Effects at Specific Values'.
pcor_si6_quad_p
sjPlot::tab_model(pcor_si6_quad)
| correct | |||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 2.49 | 2.26 – 2.75 | <0.001 |
| set size [1st degree] | 626.67 | 5.25 – 74768.23 | 0.008 |
| set size [2nd degree] | 0.00 | 0.00 – 0.00 | <0.001 |
| Random Effects | |||
| σ2 | 3.29 | ||
| τ00 ID | 0.48 | ||
| τ11 ID.poly(set_size, 2)1 | 414.12 | ||
| τ11 ID.poly(set_size, 2)2 | 366.84 | ||
| ρ01 | 0.28 | ||
| -0.48 | |||
| ICC | 0.15 | ||
| N ID | 275 | ||
| Observations | 8250 | ||
| Marginal R2 / Conditional R2 | 0.017 / 0.163 | ||
Neither linear nor quadratic signif effects
summary(pcor_p2_test_quad <-
glmer(correct ~ poly(set_size, 2) +
(poly(set_size, 2)|ID),
data=test_df %>% filter(phase==2),
family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: correct ~ poly(set_size, 2) + (poly(set_size, 2) | ID)
## Data: test_df %>% filter(phase == 2)
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 12074.6 12144.0 -6028.3 12056.6 16491
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.4912 0.1678 0.2476 0.3975 1.4058
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 1.513 1.23
## poly(set_size, 2)1 3100.386 55.68 0.02
## poly(set_size, 2)2 2531.491 50.31 0.06 -0.15
## Number of obs: 16500, groups: ID, 275
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.3021 0.0814 28.280 <2e-16 ***
## poly(set_size, 2)1 -3.1372 3.4343 -0.913 0.361
## poly(set_size, 2)2 -4.5952 3.7001 -1.242 0.214
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) p(_,2)1
## ply(st_,2)1 0.010
## ply(st_,2)2 0.025 0.069
sjPlot::tab_model(pcor_p2_test_quad)
| correct | |||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 10.00 | 8.52 – 11.72 | <0.001 |
| set size [1st degree] | 0.04 | 0.00 – 36.38 | 0.361 |
| set size [2nd degree] | 0.01 | 0.00 – 14.25 | 0.214 |
| Random Effects | |||
| σ2 | 3.29 | ||
| τ00 ID | 1.51 | ||
| τ11 ID.poly(set_size, 2)1 | 3100.39 | ||
| τ11 ID.poly(set_size, 2)2 | 2531.49 | ||
| ρ01 | 0.02 | ||
| 0.06 | |||
| ICC | 0.36 | ||
| N ID | 275 | ||
| Observations | 16500 | ||
| Marginal R2 / Conditional R2 | 0.000 / 0.361 | ||
pcor_p2_quad_p <- sjPlot::plot_model(pcor_p2_test_quad,
type = "pred", terms = c("set_size [all]")) +
ga + ap + lp + tp + xlab("Set size") +
ylab("Proportion correct") +
ggtitle("Quadratic model regression \npredictions \n Phase 2 test")
pcor_p2_quad_p
pcor_p1_test_quad_p + pcor_si6_quad_p +pcor_p2_quad_p
pcor_p2_quad_p
rl_off_p <-
ggplot(m35, aes(x=rl_off)) +
geom_vline(xintercept=median(m35$rl_off), size=4) +
geom_histogram(fill="white", color="black") + ga + ap + ylab("") +
ggtitle(TeX('$\\RL^{off}')) + tp + xlab("")
rl_off_p
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#ggsave("../paper/figs/pieces/rl_off.png", rl_off_p, height = 4, width=11, dpi=300)
little_coop <- data.frame(m35 %>% filter(rl_off < median(m35$rl_off)))$ID
high_coop <- data.frame(m35 %>% filter(rl_off > median(m35$rl_off)))$ID
p1_test_emp_trials_m35_lc <- test_df %>% filter(phase==1 & ID %in% little_coop)
pcor_p1_test_emp_m35_lc <-
data.frame(p1_test_emp_trials_m35_lc %>% group_by(set_size) %>% summarize(m=mean(correct)))
pcor_p1_test_emp_m35_lc_ID <-
data.frame(p1_test_emp_trials_m35_lc %>% group_by(set_size, ID) %>% summarize(m=mean(correct)))
## `summarise()` has grouped output by 'set_size'. You can override using the
## `.groups` argument.
p2_test_emp_trials_m35_lc <- test_df %>% filter(phase==2 & ID %in% little_coop)
pcor_p2_test_emp_m35_lc <-
data.frame(p2_test_emp_trials_m35_lc %>% group_by(set_size) %>% summarize(m=mean(correct)))
pcor_p2_test_emp_m35_lc_ID <-
data.frame(p2_test_emp_trials_m35_lc %>% group_by(set_size, ID) %>% summarize(m=mean(correct)))
## `summarise()` has grouped output by 'set_size'. You can override using the
## `.groups` argument.
si6_emp_trials_m35_lc <- learn_df %>% filter(stim_iter==6 & ID %in% little_coop)
pcor_si6_m35_lc <-
data.frame(si6_emp_trials_m35_lc %>% group_by(set_size) %>% summarize(m=mean(correct)))
pcor_si6_m35_lc_ID <-
data.frame(si6_emp_trials_m35_lc %>% group_by(set_size, ID) %>% summarize(m=mean(correct)))
## `summarise()` has grouped output by 'set_size'. You can override using the
## `.groups` argument.
emp_p1_test_m35_lc <- ggplot(pcor_p1_test_emp_m35_lc, aes(x=as.factor(set_size), y=m, fill=as.factor(set_size))) +
geom_hline(yintercept = seq(.1, 1, .1), alpha=.3) +
ga + ap + xlab("Set size") + ylab("Proportion correct") + tol + ylim(-.022, 1.022) + tp +
geom_jitter(data=pcor_p1_test_emp_m35_lc_ID, width=.17, height=.02, alpha=.8, aes(color=as.factor(set_size))) +
geom_bar(stat="identity", color="black", alpha=.8) +
tp + ggtitle("Test phase 1")
emp_p2_test_m35_lc <- ggplot(pcor_p2_test_emp_m35_lc, aes(x=as.factor(set_size), y=m, fill=as.factor(set_size))) +
geom_hline(yintercept = seq(.1, 1, .1), alpha=.3) +
ga + ap + xlab("Set size") + ylab("Proportion correct") + tol + ylim(-.022, 1.022) + tp +
geom_jitter(data=pcor_p2_test_emp_m35_lc_ID, width=.17, height=.02, alpha=.8, aes(color=as.factor(set_size))) +
geom_bar(stat="identity", color="black", alpha=.8) +
tp + ggtitle("Test phase 2")
emp_si6_lc <-
ggplot(pcor_si6_m35_lc, aes(x=as.factor(set_size), y=m, fill=as.factor(set_size))) +
geom_hline(yintercept = seq(.1, 1, .1), alpha=.3) +
ga + ap + xlab("Set size") + ylab("Proportion correct") + tol + ylim(-.022, 1.022) + tp +
geom_jitter(data=pcor_si6_m35_lc_ID, width=.17, height=.02, alpha=.8, aes(color=as.factor(set_size))) +
geom_bar(stat="identity", color="black", alpha=.8) +
tp + ggtitle("Stimulus iteration 6")
low_blunt <- emp_p1_test_m35_lc + emp_si6_lc + emp_p2_test_m35_lc
High blunters
p1_test_emp_trials_m35_hc <- test_df %>% filter(phase==1 & ID %in% high_coop)
pcor_p1_test_emp_m35_hc <-
data.frame(p1_test_emp_trials_m35_hc %>% group_by(set_size) %>% summarize(m=mean(correct)))
pcor_p1_test_emp_m35_hc_ID <-
data.frame(p1_test_emp_trials_m35_hc %>% group_by(set_size, ID) %>% summarize(m=mean(correct)))
## `summarise()` has grouped output by 'set_size'. You can override using the
## `.groups` argument.
p2_test_emp_trials_m35_hc <- test_df %>% filter(phase==2 & ID %in% high_coop)
pcor_p2_test_emp_m35_hc <-
data.frame(p2_test_emp_trials_m35_hc %>% group_by(set_size) %>% summarize(m=mean(correct)))
pcor_p2_test_emp_m35_hc_ID <-
data.frame(p2_test_emp_trials_m35_hc %>% group_by(set_size, ID) %>% summarize(m=mean(correct)))
## `summarise()` has grouped output by 'set_size'. You can override using the
## `.groups` argument.
si6_emp_trials_m35_hc <- learn_df %>% filter(stim_iter==6 & ID %in% high_coop)
pcor_si6_m35_hc <-
data.frame(si6_emp_trials_m35_hc %>% group_by(set_size) %>% summarize(m=mean(correct)))
pcor_si6_m35_hc_ID <-
data.frame(si6_emp_trials_m35_hc %>% group_by(set_size, ID) %>% summarize(m=mean(correct)))
## `summarise()` has grouped output by 'set_size'. You can override using the
## `.groups` argument.
emp_p1_test_m35_hc <- ggplot(pcor_p1_test_emp_m35_hc, aes(x=as.factor(set_size), y=m, fill=as.factor(set_size))) +
geom_hline(yintercept = seq(.1, 1, .1), alpha=.3) +
ga + ap + xlab("Set size") + ylab("Proportion correct") + tol + ylim(-.022, 1.022) + tp +
geom_jitter(data=pcor_p1_test_emp_m35_hc_ID, width=.17, height=.02, alpha=.8, aes(color=as.factor(set_size))) +
geom_bar(stat="identity", color="black", alpha=.8) #+
#tp + ggtitle("Phase 1")
emp_p2_test_m35_hc <- ggplot(pcor_p2_test_emp_m35_hc, aes(x=as.factor(set_size), y=m, fill=as.factor(set_size))) +
geom_hline(yintercept = seq(.1, 1, .1), alpha=.3) +
ga + ap + xlab("Set size") + ylab("Proportion correct") + tol + ylim(-.022, 1.022) + tp +
geom_jitter(data=pcor_p2_test_emp_m35_hc_ID, width=.17, height=.02, alpha=.8, aes(color=as.factor(set_size))) +
geom_bar(stat="identity", color="black", alpha=.8) #+
# tp + ggtitle("Phase 2")
emp_si6_hc <-
ggplot(pcor_si6_m35_hc, aes(x=as.factor(set_size), y=m, fill=as.factor(set_size))) +
geom_hline(yintercept = seq(.1, 1, .1), alpha=.3) +
ga + ap + xlab("Set size") + ylab("Proportion correct") + tol + ylim(-.022, 1.022) + tp +
geom_jitter(data=pcor_si6_m35_hc_ID, width=.17, height=.02, alpha=.8, aes(color=as.factor(set_size))) +
geom_bar(stat="identity", color="black", alpha=.8) #+
#tp + ggtitle("Stimulus iteration 6")
high_blunt <- emp_p1_test_m35_hc + emp_si6_hc + emp_p2_test_m35_hc
low_vs_high <- low_blunt/high_blunt
low_vs_high
#ggsave("../paper/figs/pieces/fig4_lohiblunt.png", low_vs_high, height = 9, width=11, dpi=300)
Put RL off into the test df
unique_ids <- unique(test_df$ID)
for (i in 1:length(unique_ids)) {
test_df[test_df$ID==unique_ids[i], "rl_off"] <- m35[m35$ID==unique_ids[i], "rl_off"]
}
# Spot check
# test_df %>% filter(ID == 24) %>% select(rl_off)
# m35 %>% filter(ID == 24) %>% select(rl_off)
# test_df %>% filter(ID == 85) %>% select(rl_off)
# m35 %>% filter(ID == 85) %>% select(rl_off)
# test_df %>% filter(ID == 2) %>% select(rl_off)
# m35 %>% filter(ID == 2) %>% select(rl_off)
# test_df %>% filter(ID == 219) %>% select(rl_off)
# m35 %>% filter(ID == 219) %>% select(rl_off)
m_rlo <- mean(m35$rl_off)
sd_up <- m_rlo +sd(m35$rl_off) # (.5*sd(m35$rl_off))
sd_down <- m_rlo - sd(m35$rl_off)# (.5*sd(m35$rl_off))
cat("\nMean", m_rlo, "\nOne SD up", sd_up, "\nOne SD down", sd_down)
##
## Mean 0.544767
## One SD up 0.8678609
## One SD down 0.221673
#quantile(m35$rl_off, seq(1/3, 1, 1/3))
#ggsave("../paper/figs/pieces/fig4_mod_rl_off_p.png", mod_rl_off_p, height = 4, width=8, dpi=300)
summary(pcor_p2_moderation <-
glmer(correct ~ scale(set_size)*scale(rl_off) +
(scale(set_size)*scale(rl_off)|ID),
data=test_df %>% filter(phase==2),
family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: correct ~ scale(set_size) * scale(rl_off) + (scale(set_size) *
## scale(rl_off) | ID)
## Data: test_df %>% filter(phase == 2)
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 12009.4 12117.4 -5990.7 11981.4 16486
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -8.6425 0.1503 0.2565 0.4166 1.5443
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 1.20709 1.0987
## scale(set_size) 0.04156 0.2039 0.19
## scale(rl_off) 0.25019 0.5002 -0.34 0.41
## scale(set_size):scale(rl_off) 0.06133 0.2476 0.68 0.41 0.41
## Number of obs: 16500, groups: ID, 275
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.32737 0.08506 27.362 < 2e-16 ***
## scale(set_size) -0.03490 0.04117 -0.848 0.397
## scale(rl_off) -0.54108 0.09283 -5.829 5.59e-09 ***
## scale(set_size):scale(rl_off) 0.45056 0.04696 9.594 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) scl(s_) scl(r_)
## scal(st_sz) -0.008
## scal(rl_ff) -0.306 0.354
## scl(s_):(_) 0.352 -0.404 -0.038
car::vif(pcor_p2_moderation)
## scale(set_size) scale(rl_off)
## 1.384028 1.160176
## scale(set_size):scale(rl_off)
## 1.212803
mod_rl_off_p <- sjPlot::plot_model(pcor_p2_moderation, type = "pred",
terms = c("set_size [all]", "rl_off [.2217, .5447, .8678]")) +
ga + ap + lp + tp + xlab("Set size") +
ylab("") +
ggtitle("Moderation of set size effect in final test phase \n by RL blunting") +
scale_color_manual(values=c("blue", "gray40", "red"), labels=c("Low", "Medium", "High")) +
scale_fill_manual(values=c("blue", "gray40", "red"), labels=c("Low", "Medium", "High"))
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
mod_rl_off_p
sjPlot::tab_model(pcor_p2_moderation)
| correct | |||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 10.25 | 8.68 – 12.11 | <0.001 |
| set size | 0.97 | 0.89 – 1.05 | 0.397 |
| rl off | 0.58 | 0.49 – 0.70 | <0.001 |
| set size × rl off | 1.57 | 1.43 – 1.72 | <0.001 |
| Random Effects | |||
| σ2 | 3.29 | ||
| τ00 ID | 1.21 | ||
| τ11 ID.scale(set_size) | 0.04 | ||
| τ11 ID.scale(rl_off) | 0.25 | ||
| τ11 ID.scale(set_size):scale(rl_off) | 0.06 | ||
| ρ01 | 0.19 | ||
| -0.34 | |||
| 0.68 | |||
| ICC | 0.32 | ||
| N ID | 275 | ||
| Observations | 16500 | ||
| Marginal R2 / Conditional R2 | 0.093 / 0.385 | ||
Low vs. high low set sizes percents reported in paper
mean(pcor_p1_test_emp_m35_hc[1:3, "m"])
## [1] 0.5612328
mean(pcor_p1_test_emp_m35_lc[1:3, "m"])
## [1] 0.729927
mean(pcor_si6_m35_hc[1:3, "m"])
## [1] 0.5788727
mean(pcor_si6_m35_lc[1:3, "m"])
## [1] 0.7226277
mean(pcor_p2_test_emp_m35_hc[1:3, "m"])
## [1] 0.7650041
mean(pcor_p2_test_emp_m35_lc[1:3, "m"])
## [1] 0.915146
Low vs. high high set sizes
mean(pcor_p1_test_emp_m35_hc[4:5, "m"])
## [1] 0.6972628
mean(pcor_p1_test_emp_m35_lc[4:5, "m"])
## [1] 0.7048814
mean(pcor_si6_m35_hc[4:5, "m"])
## [1] 0.6803832
mean(pcor_si6_m35_lc[4:5, "m"])
## [1] 0.7076642
mean(pcor_p2_test_emp_m35_hc[4:5, "m"])
## [1] 0.8436131
mean(pcor_p2_test_emp_m35_lc[4:5, "m"])
## [1] 0.8611314
Checking low vs. high asymptotic performance — set sizes 1 and 2 are very comparable
lc_asym <- learn_df %>% filter(stim_iter==10 & ID %in% little_coop) %>% group_by(set_size) %>% summarize(m=mean(correct))
hc_asym <- learn_df %>% filter(stim_iter==10 & ID %in% high_coop) %>% group_by(set_size) %>% summarize(m=mean(correct))
mean(unlist(lc_asym[1:3, "m"]))
## [1] 0.9578264
mean(unlist(hc_asym[1:3, "m"]))
## [1] 0.9551906
mean(lc_asym[1:3, "m"])
## Warning in mean.default(lc_asym[1:3, "m"]): argument is not numeric or logical:
## returning NA
## [1] NA
RL off is not much correlated with WM pars
ComparePars(m35$rl_off, m35$phi, use_identity_line = 0)
ComparePars(m35$rl_off, m35$kappa, use_identity_line = 0)
ComparePars(m35$rl_off, m35$rho, use_identity_line = 0)
A key change in this task from the classic RLWM was having punishment | neutral | reward bandit arms instead of just reward | neutral | neutral. But it’s possible that pts would focus on maximizing reward and not make much distinguish punishment vs. neutral. The below examine neutral > punishment preference during learning and test (retention).
learn_error_df <- learn_df %>% filter(correct==0)
np_summs_m <- learn_error_df %>% group_by(stim_iter) %>%
summarize(mw=mean(worst), mn=mean(neutral), n=n()) %>% mutate(neutral_pref=(mn-mw))
# Within-subject adjusted mean and raw don't agree presumably because of amount of
# variation given sparsity of errors/subject so will just use the raw means for plotting
# so can do accurate compare to model preds
# np_si <- Rmisc::summarySEwithin(np_summs,
# measurevar = "neutral_pref",
# withinvars = c("stim_iter"),
# idvar = "ID")
# These do disagree before this decimal point
# round(np_summs_m$neutral_pref, 4)
# round(np_si$neutral_pref, 4)
emp_np_si_plot <-
ggplot(np_summs_m, aes(x=stim_iter, y=neutral_pref, fill=as.numeric(stim_iter))) +
geom_hline(yintercept=.0, size=1.5, color="gray57") + # chance line
geom_bar(stat="identity", color="black") +
#geom_errorbar(aes(ymin=neutral_pref-se, ymax=neutral_pref+se), width=.2) +
geom_bar(data=np_summs_m,
aes(x=stim_iter, y=neutral_pref, fill=as.numeric(stim_iter)),
stat="identity", color="black", alpha=.85) +
annotate("rect", xmin=5.5, xmax=10.5, ymin=0, ymax=.22, alpha=0.2, fill="gray57") +
annotate("text", x=3, y=.20, label="Phase 1", size=8) +
annotate("text", x=7.5, y=.20, label="Phase 2", size=8) +
ga + ap + tol + xlab("Simulus iteration") + ylab("") + tp +
ggtitle("Empirical") + scale_color_gradient2() + ylim(-.1, .23)
emp_np_si_plot
#move_layers(emp_np_si_plot, "GeomPoint", position = "top")
Capture by sims
np_si_sims_summs <- m35_s %>%
filter(type=="learning" & corrects==0) %>% group_by(stim_iter) %>% #filter(stim_iter %in% c(2:10)) %>%
summarize(mw=mean(worsts), mn=mean(neutrals), n=n()) %>% mutate(neutral_pref=(mn-mw))
np_si_sims_summs_var <- m35_s %>%
filter(type=="learning" & corrects==0) %>% group_by(stim_iter, iter) %>% #filter(stim_iter %in% c(2:10)) %>%
summarize(mw=mean(worsts), mn=mean(neutrals), n=n()) %>% mutate(neutral_pref=(mn-mw))
## `summarise()` has grouped output by 'stim_iter'. You can override using the
## `.groups` argument.
m35_np_si_plot <- ggplot(np_si_sims_summs,
aes(x=stim_iter, y=neutral_pref, fill=as.numeric(stim_iter))) +
geom_hline(yintercept=.0, size=1.5, color="gray57") + # chance line
geom_bar(stat="identity", color="black") +
geom_jitter(data=np_si_sims_summs_var, aes(x=stim_iter, y=neutral_pref), height=0,
width=.2, pch=21) +
annotate("rect", xmin=5.5, xmax=10.5, ymin=0, ymax=.22, alpha=0.2, pch=21) +
annotate("text", x=3, y=.20, label="Phase 1", size=8) +
annotate("text", x=7.5, y=.20, label="Phase 2", size=8) +
ga + ap + tol + xlab("Simulus iteration") + ylab("") + tp +
ylim(-.05, .23) +
ggtitle("Simulated") + scale_color_gradient2() +
scale_x_continuous(breaks=seq(1, 10, 1)) + ylim(-.1, .23)
## Warning in annotate("rect", xmin = 5.5, xmax = 10.5, ymin = 0, ymax = 0.22, :
## Ignoring unknown parameters: `shape`
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
m35_np_si_plot
Simulation w no punish bonus off
m35_luu <- read.csv("../model_res/sims/SIM_RunRLWMPRewLesionNotPun67064.csv") # Updated non-bug 8-1-24 model
np_si_sims_summs_luu <- m35_luu %>%
filter(type=="learning" & corrects==0) %>% group_by(stim_iter) %>%
summarize(mw=mean(worsts), mn=mean(neutrals), n=n()) %>% mutate(neutral_pref=(mn-mw))
np_si_sims_summs_var_luu <- m35_luu %>%
filter(type=="learning" & corrects==0) %>% group_by(stim_iter, iter) %>%
summarize(mw=mean(worsts), mn=mean(neutrals), n=n()) %>% mutate(neutral_pref=(mn-mw))
## `summarise()` has grouped output by 'stim_iter'. You can override using the
## `.groups` argument.
m35_np_si_luu_plot <-
ggplot(np_si_sims_summs_luu,
aes(x=stim_iter, y=neutral_pref, fill=as.numeric(stim_iter))) +
geom_hline(yintercept=.0, size=1.5, color="gray57") + # chance line
geom_bar(stat="identity", color="black") +
geom_jitter(data=np_si_sims_summs_var_luu, aes(x=stim_iter, y=neutral_pref), pch=21, width=.2) +
annotate("rect", xmin=5.5, xmax=10.5, ymin=0, ymax=.22, alpha=0.2, fill="gray57") +
annotate("text", x=3, y=.20, label="Phase 1", size=8) +
annotate("text", x=7.5, y=.20, label="Phase 2", size=8) +
ga + ap + tol + xlab("Simulus iteration") + ylab("") + tp + ylim(-.1, .23) +
ggtitle(TeX("Sim.: Non-pun. bonus = 0")) + scale_color_gradient2() +
scale_x_continuous(breaks=seq(1, 10, 1)) + theme(plot.title = element_text(size = 18))
#theme(plot.title = element_text(vjust = -2.5))
m35_np_si_luu_plot
np_si <- emp_np_si_plot + m35_np_si_plot + m35_np_si_luu_plot
np_si
#ggsave("../paper/figs/pieces/fig5_neut-pref_emp-and-sims.png", np_si, height = 6, width=12, dpi=300)
# Didn't converge, so do Bayesian version
# summary(pneut_phase1_le <- glmer(neutral ~ scale(stim_iter) + (0 + scale(stim_iter) |ID),
# data=learn_error_df %>% filter(phase==1), family="binomial", control = glmerControl(optimizer = "bobyqa")))
# Now runnning in brms_s script
# pneut_time <- brm(
# neutral ~ scale(stim_iter) + (1 |ID),
# data = learn_error_df %>% filter(phase==1),
# family = bernoulli(link = "logit"),
# warmup=2e3,
# iter=4e3,
# chains=5,
# cores=5,
# control= list(adapt_delta = 0.9))
Note: brms files not shared in public-data due to file
size but can be reproduced using brms-s.R
np_p1_brms <-
read.csv("../model_res/brms_res/neutral_pref_time_effect_phase1__66703__.csv")
np_p2_brms <-
read.csv("../model_res/brms_res/neutral_pref_time_effect_phase2__28768__.csv")
hist(np_p1_brms$b_scalestim_iter, breaks=100)
hist(np_p2_brms$b_scalestim_iter, breaks=100)
test_error_df <- test_df %>% filter(correct==0)
# np_test_summs <- test_error_df %>% group_by(ID, phase) %>%
# summarize(mw=mean(worst), mn=mean(neutral), n=n()) %>% mutate(neutral_pref=(mn-mw))
np_test_summs_m <- test_error_df %>% group_by(phase) %>%
summarize(mw=mean(worst), mn=mean(neutral), n=n()) %>% mutate(neutral_pref=(mn-mw))
# np_test_si <- Rmisc::summarySEwithin(np_test_summs,
# measurevar = "neutral_pref",
# withinvars = c("phase"),
# idvar = "ID")
#
# # Within-subject adjusted mean and raw don't agree presumably because of amount of
# # variation given sparsity of errors/subject so will just use the raw means for plotting
# # so can do accurate compare to model preds
# np_test_si_sew <- Rmisc::summarySEwithin(np_test_summs,
# measurevar = "neutral_pref",
# withinvars = c("phase"),
# idvar = "ID")
#
# round(np_test_summs_m$neutral_pref, 4)
# round(np_test_si_sew$neutral_pref, 4)
In contrast to during learning, participants do not appear to retain much ability to avoid the wost stimuli during test
np_test_summs_m
## # A tibble: 2 × 5
## phase mw mn n neutral_pref
## <int> <dbl> <dbl> <int> <dbl>
## 1 1 0.494 0.506 5129 0.0119
## 2 2 0.488 0.512 2458 0.0236
emp_np_si_test_plot <-
ggplot(np_test_summs_m, aes(x=as.factor(phase), y=neutral_pref, fill=as.factor(phase))) +
geom_hline(yintercept=.0, size=1.5, color="gray57") + # chance line
geom_bar(stat="identity", color="black") +
#geom_errorbar(aes(ymin=neutral_pref-se, ymax=neutral_pref+se), width=.2) +
ga + ap + tol + xlab("Phase") + ylab("") + tp +
ggtitle(#"Neutral preference during test
"Empirical") + scale_fill_manual(values=c("gray81", "gray40")) + ylim(-.05, .23)
emp_np_si_test_plot
# ggsave("../paper/figs/pieces/fig6_emp_np_si.png", emp_np_si_test_plot, height = 4, width=5, dpi=300)
And statistically there is no significant ability to pick neutral over worst, or better ability to do so by the second versus first phase.
(These are not singular issue so didn’t need to fit Bayesian model)
summary(pneut_test <- glmer(neutral ~ 1 + (1|ID),
data=test_error_df, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: neutral ~ 1 + (1 | ID)
## Data: test_error_df
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 10474.7 10488.6 -5235.4 10470.7 7585
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.4510 -0.9779 0.7725 0.9626 1.3129
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0.09459 0.3076
## Number of obs: 7587, groups: ID, 275
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.03924 0.03108 1.263 0.207
summary(pneut_test_phase <- glmer(neutral ~ phase + (phase|ID),
data=test_error_df, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: neutral ~ phase + (phase | ID)
## Data: test_error_df
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 10477.0 10511.6 -5233.5 10467.0 7582
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.6292 -0.9788 0.7230 0.9641 1.3537
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 0.12691 0.3562
## phase 0.06782 0.2604 -0.61
## Number of obs: 7587, groups: ID, 275
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.02154 0.07496 -0.287 0.774
## phase 0.04757 0.05488 0.867 0.386
##
## Correlation of Fixed Effects:
## (Intr)
## phase -0.910
Bayesian version for more direct comparison to learning where only Bayesian model converged
np_test_brms <- read.csv("../model_res/brms_res/neutral_pref_test_effect__33330__.csv")
hist(np_test_brms$b_Intercept, breaks=100)
However store the frequentist REs for correlating with test diffs
test_res <- data.frame(ranef(pneut_test)$ID)
The lack of retained ability to avoid punishment over neutral is unsurprising because the RL learning rate from negative PEs fit much lower than that from positive PEs
cat("\nAlpha for pos. PEs quantiles: \n")
##
## Alpha for pos. PEs quantiles:
round(quantile(m35$alpha_pos), 6)
## 0% 25% 50% 75% 100%
## 0.000767 0.005559 0.007560 0.010339 0.048076
cat("\nAlpha for neg. PEs quantiles: \n")
##
## Alpha for neg. PEs quantiles:
round(quantile(m35$alpha_neg), 6)
## 0% 25% 50% 75% 100%
## 0.000000 0.000000 0.000122 0.002330 0.023870
Take log given strong negative skew
hist(m35$alpha_neg, breaks=100)
hist(log(m35$alpha_neg), breaks=100)
hist(m35$alpha_pos, breaks=100)
hist(log(m35$alpha_pos), breaks=100)
t.test(log(m35$alpha_neg), log(m35$alpha_pos), paired = TRUE)
##
## Paired t-test
##
## data: log(m35$alpha_neg) and log(m35$alpha_pos)
## t = -16.124, df = 274, p-value < 2.2e-16
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
## -7.853473 -6.144423
## sample estimates:
## mean difference
## -6.998948
m35_test_sim_error_df <- m35_s_test %>% filter(corrects==0)
np_test_summs_m35_sim <- m35_test_sim_error_df %>% group_by(phase) %>%
summarize(mw=mean(worsts), mn=mean(neutrals), n=n()) %>% mutate(neutral_pref=(mn-mw))
np_test_summs_m35_sim_var <- m35_test_sim_error_df %>% group_by(phase, iter) %>%
summarize(mw=mean(worsts), mn=mean(neutrals), n=n()) %>% mutate(neutral_pref=(mn-mw))
## `summarise()` has grouped output by 'phase'. You can override using the
## `.groups` argument.
And the model correctly predicts poor retained ability to avoid punishment over neutral overall
sim_m35_np_si_test_plot <-
ggplot(np_test_summs_m35_sim , aes(x=as.factor(phase), y=neutral_pref, fill=as.factor(phase))) +
geom_hline(yintercept=.0, size=1.5, color="gray57") + # chance line
geom_jitter(data=np_test_summs_m35_sim_var, height=0, width=.2, size=3, alpha=.8) +
geom_bar(stat="identity", color="black", alpha=.8) +
ga + ap + tol + xlab("Phase") + ylab("") + tp +
ggtitle(#"Neutral preference during Test
"Simulated") + scale_fill_manual(values=c("gray81", "gray40")) + ylim(-.05, .23)
sim_m35_np_si_test_plot
Whereas a model that sets \(\alpha_{neg} = \alpha_{pos}\) incorrectly predicts a retained neutral preference at a similar level to what is observed empirically at the highest points during learning
m35_s_lesionalphaneg <-
read.csv("../model_res/sims/SIM_RunRLWMPRewLesionAlphaNeg27102.csv")
m35_test_sim_error_df_lan <- m35_s_lesionalphaneg %>% filter(corrects==0)
np_test_summs_m35_sim_lan <- m35_test_sim_error_df_lan %>% group_by(phase) %>%
summarize(mw=mean(worsts), mn=mean(neutrals), n=n()) %>% mutate(neutral_pref=(mn-mw))
np_test_summs_m35_sim_var_lan <- m35_test_sim_error_df_lan %>% group_by(phase, iter) %>%
summarize(mw=mean(worsts), mn=mean(neutrals), n=n()) %>% mutate(neutral_pref=(mn-mw))
## `summarise()` has grouped output by 'phase'. You can override using the
## `.groups` argument.
sim_m35_np_si_test_plot_lan <-
ggplot(np_test_summs_m35_sim_lan , aes(x=as.factor(phase), y=neutral_pref, fill=as.factor(phase))) +
geom_hline(yintercept=.0, size=1.5, color="gray57") + # chance line
geom_jitter(data=np_test_summs_m35_sim_var_lan, height=0, width=.2, size=3, alpha=.8) +
geom_bar(stat="identity", color="black", alpha=.8) +
ga + ap + tol + xlab("Phase") + ylab("") + tp +
ggtitle(TeX("Simulated: $\\alpha^{-} = \\alpha^{+}$")) +
scale_fill_manual(values=c("gray81", "gray40")) + ylim(-.05, .23)
#ggtitle("Neutral preference during Test \nSimulated equating TeX('$\\alpha^{-}')") +
sim_m35_np_si_test_plot_lan
test_np_plots <-
emp_np_si_test_plot + sim_m35_np_si_test_plot #+ sim_m35_np_si_test_plot_lan
test_np_plots
ggsave("../paper/figs/pieces/fig6_TEST_neut-pref_emp-and-sims.png", test_np_plots, height = 6, width=12, dpi=300)
Code if alpha pos is greater
m35$alpha_pos_greater <- if_else(m35$alpha_pos > m35$alpha_neg, 1, 0)
table(m35$alpha_pos_greater)
##
## 0 1
## 17 258
table(m35$alpha_pos_greater)[2]/sum(table(m35$alpha_pos_greater))
## 1
## 0.9381818
alpha_comp <- ggplot(m35,
aes(x=alpha_neg, alpha_pos, fill=as.factor(alpha_pos_greater))) +
geom_point(size=4, pch=21) +
ga + ap + tp +
stp + tol + ylim(0, .05) + xlim(0, .05) +
geom_line(aes(y=alpha_neg)) +
xlab(TeX("$\\alpha^{-} ")) +
ylab(TeX("$\\alpha^{+} ")) + scale_fill_manual(values=c("purple", "orange")) +
theme(axis.title = element_text(size=50))
# ggtitle(model_char, subtitle=str) +
# ylab(ychar) + xlab(xchar)
alpha_comp
alpha_and_sim <- alpha_comp + sim_m35_np_si_test_plot_lan
# ggsave("../paper/figs/pieces/fig5_alpha_and_sim.png", alpha_and_sim, height = 6, width=12, dpi=300)
However, there are some individual differences — with those fitting a higher learning from negative PEs showing better ability to retain a neutral preference (and note this correlation is expected to be attenuated due to measurement error due to the relatively poor parameter recovery for the learning rate for negative PEs shown earlier)
hist(m35$alpha_neg)
ComparePars(m35$alpha_neg, test_res$X.Intercept.,
model_char="", xchar = "Learning rate from negative RPEs",
ychar = "Avoidance of punishment \n over neutral at test",
use_identity_line = 0)
ComparePars(log(m35$alpha_neg), test_res$X.Intercept.,
model_char="", xchar = "Log (learning rate from negative RPEs)",
ychar = "Avoidance of punishment \n over neutral at test",
use_identity_line = 0)
# qdf %>% filter(catch_q_1 != 1) %>% select(ID)
# qdf %>% filter(catch_q_2 != 1) %>% select(ID)
qdf <- read.csv("../data/questionnaire_df.csv")
demogs <- read.csv("../data/demogs_deident.csv")
if (all(demogs$ID==qdf$ID)) qdf$group <- demogs$group
qdf <- qdf %>% filter(!ID==255)
demogs <- demogs %>% filter(!ID==255)
learn_df <- learn_df %>% filter(!ID==255)
test_df <- test_df %>% filter(!ID==255)
Tiny amount of missing cases so fine to mean impute
length(which(is.na(qdf %>% select(contains("BDI"))))) # 3 missing BDI
## [1] 3
length(which(is.na(qdf %>% select(contains("GAD"))))) # 4 GAD
## [1] 4
length(which(is.na(qdf %>% select(contains("RRS"))))) # 3 RRS
## [1] 3
qdf_before_imp <- qdf
for (i in 1:nrow(qdf)) {
if (any(is.na(qdf[i, grep("BDI", names(qdf))]))) {
this_row <- qdf[i, grep("BDI", names(qdf))]
na_idx <- which(is.na(qdf[i, grep("BDI", names(qdf))]))
this_row[na_idx] <- mean(as.numeric(this_row), na.rm=TRUE)
qdf[i, grep("BDI", names(qdf))] <- this_row
}
}
for (i in 1:nrow(qdf)) {
if (any(is.na(qdf[i, grep("GAD", names(qdf))]))) {
this_row <- qdf[i, grep("GAD", names(qdf))]
na_idx <- which(is.na(qdf[i, grep("GAD", names(qdf))]))
this_row[na_idx] <- mean(as.numeric(this_row), na.rm=TRUE)
qdf[i, grep("GAD", names(qdf))] <- this_row
}
}
for (i in 1:nrow(qdf)) {
if (any(is.na(qdf[i, grep("RRS", names(qdf))]))) {
this_row <- qdf[i, grep("RRS", names(qdf))]
na_idx <- which(is.na(qdf[i, grep("RRS", names(qdf))]))
this_row[na_idx] <- mean(as.numeric(this_row), na.rm=TRUE)
qdf[i, grep("RRS", names(qdf))] <- this_row
}
}
Given tiny amount of missing cases so shouldn’t skew values much, doing NA rm for pres
length(which(is.na(qdf %>% select(contains("BDI"))))) # 3 missing BDI
## [1] 0
length(which(is.na(qdf %>% select(contains("GAD"))))) # 4 GAD
## [1] 0
length(which(is.na(qdf %>% select(contains("RRS"))))) # 3 RRS
## [1] 0
Sanity check unchanged
# ComparePars(rowSums(qdf[, grep("RRS", names(qdf))]),
# rowSums(qdf_before_imp[, grep("RRS", names(qdf_before_imp))]))
# ComparePars(rowSums(qdf[, grep("BDI", names(qdf))]),
# rowSums(qdf_before_imp[, grep("BDI", names(qdf_before_imp))]))
# ComparePars(rowSums(qdf[, grep("GAD", names(qdf))]),
# rowSums(qdf_before_imp[, grep("GAD", names(qdf_before_imp))]))
Reliability stats
psych::alpha(qdf[, grep("BDI", names(qdf))])
##
## Reliability analysis
## Call: psych::alpha(x = qdf[, grep("BDI", names(qdf))])
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.96 0.96 0.97 0.54 25 0.0033 0.68 0.66 0.54
##
## 95% confidence boundaries
## lower alpha upper
## Feldt 0.95 0.96 0.97
## Duhachek 0.95 0.96 0.97
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## BDI.1 0.96 0.96 0.97 0.54 24 0.0034 0.0107 0.54
## BDI.2 0.96 0.96 0.97 0.54 23 0.0035 0.0102 0.54
## BDI.3 0.96 0.96 0.97 0.54 23 0.0036 0.0099 0.54
## BDI.4 0.96 0.96 0.97 0.54 23 0.0035 0.0102 0.54
## BDI.5 0.96 0.96 0.97 0.55 24 0.0034 0.0111 0.55
## BDI.6 0.96 0.96 0.97 0.55 24 0.0034 0.0107 0.55
## BDI.7 0.96 0.96 0.97 0.54 23 0.0036 0.0099 0.54
## BDI.8 0.96 0.96 0.97 0.54 23 0.0036 0.0101 0.54
## BDI.9 0.96 0.96 0.97 0.55 24 0.0034 0.0103 0.54
## BDI.10 0.96 0.96 0.97 0.55 24 0.0034 0.0107 0.56
## BDI.11 0.96 0.96 0.97 0.54 24 0.0034 0.0111 0.55
## BDI.12 0.96 0.96 0.97 0.54 23 0.0036 0.0100 0.54
## BDI.13 0.96 0.96 0.97 0.54 24 0.0035 0.0109 0.54
## BDI.14 0.96 0.96 0.97 0.54 23 0.0036 0.0098 0.54
## BDI.15 0.96 0.96 0.97 0.54 23 0.0035 0.0104 0.54
## BDI.16 0.96 0.96 0.97 0.55 25 0.0033 0.0099 0.56
## BDI.17 0.96 0.96 0.97 0.55 24 0.0034 0.0106 0.55
## BDI.18 0.96 0.96 0.97 0.55 25 0.0033 0.0099 0.56
## BDI.19 0.96 0.96 0.97 0.54 24 0.0034 0.0106 0.54
## BDI.20 0.96 0.96 0.97 0.54 23 0.0035 0.0107 0.54
## BDI.21 0.96 0.96 0.97 0.55 25 0.0033 0.0096 0.56
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## BDI.1 275 0.73 0.74 0.72 0.71 0.56 0.80
## BDI.2 275 0.80 0.79 0.78 0.77 0.85 0.95
## BDI.3 275 0.84 0.83 0.83 0.81 0.93 1.03
## BDI.4 275 0.84 0.84 0.83 0.82 0.71 0.84
## BDI.5 275 0.73 0.73 0.71 0.70 0.60 0.77
## BDI.6 275 0.72 0.72 0.70 0.69 0.56 0.95
## BDI.7 275 0.84 0.83 0.83 0.81 0.88 1.07
## BDI.8 275 0.84 0.83 0.83 0.81 0.76 0.94
## BDI.9 275 0.70 0.70 0.69 0.67 0.37 0.66
## BDI.10 275 0.68 0.68 0.65 0.64 0.47 0.89
## BDI.11 275 0.73 0.74 0.73 0.71 0.56 0.74
## BDI.12 275 0.85 0.85 0.85 0.83 0.78 0.90
## BDI.13 275 0.78 0.78 0.77 0.76 0.52 0.85
## BDI.14 275 0.83 0.83 0.83 0.81 0.75 1.02
## BDI.15 275 0.80 0.80 0.79 0.77 0.89 0.92
## BDI.16 275 0.62 0.62 0.59 0.58 0.88 0.94
## BDI.17 275 0.70 0.71 0.69 0.67 0.51 0.75
## BDI.18 275 0.61 0.62 0.59 0.57 0.53 0.80
## BDI.19 275 0.74 0.74 0.73 0.71 0.65 0.84
## BDI.20 275 0.81 0.81 0.81 0.79 0.79 0.90
## BDI.21 275 0.61 0.61 0.58 0.57 0.64 0.93
##
## Non missing response frequency for each item
## 0 0.5 1 1.05 1.35 2 3 miss
## BDI.1 0.58 0 0.32 0 0 0.05 0.05 0
## BDI.2 0.45 0 0.34 0 0 0.13 0.08 0
## BDI.3 0.46 0 0.26 0 0 0.17 0.11 0
## BDI.4 0.51 0 0.30 0 0 0.16 0.03 0
## BDI.5 0.56 0 0.32 0 0 0.10 0.03 0
## BDI.6 0.68 0 0.17 0 0 0.07 0.09 0
## BDI.7 0.52 0 0.21 0 0 0.16 0.12 0
## BDI.8 0.53 0 0.24 0 0 0.17 0.06 0
## BDI.9 0.71 0 0.24 0 0 0.03 0.02 0
## BDI.10 0.72 0 0.17 0 0 0.03 0.08 0
## BDI.11 0.57 0 0.32 0 0 0.09 0.02 0
## BDI.12 0.48 0 0.30 0 0 0.16 0.05 0
## BDI.13 0.65 0 0.23 0 0 0.06 0.06 0
## BDI.14 0.58 0 0.19 0 0 0.13 0.10 0
## BDI.15 0.42 0 0.32 0 0 0.19 0.06 0
## BDI.16 0.43 0 0.32 0 0 0.17 0.07 0
## BDI.17 0.61 0 0.28 0 0 0.08 0.03 0
## BDI.18 0.63 0 0.25 0 0 0.08 0.04 0
## BDI.19 0.56 0 0.27 0 0 0.14 0.03 0
## BDI.20 0.47 0 0.33 0 0 0.13 0.06 0
## BDI.21 0.61 0 0.20 0 0 0.12 0.07 0
psych::alpha(qdf[, grep("GAD", names(qdf))])
##
## Reliability analysis
## Call: psych::alpha(x = qdf[, grep("GAD", names(qdf))])
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.94 0.94 0.94 0.67 15 0.0058 0.84 0.83 0.66
##
## 95% confidence boundaries
## lower alpha upper
## Feldt 0.92 0.94 0.95
## Duhachek 0.92 0.94 0.95
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## GAD_1 0.92 0.92 0.92 0.66 11 0.0074 0.0074 0.65
## GAD_2 0.92 0.92 0.91 0.66 11 0.0074 0.0059 0.65
## GAD_3 0.92 0.92 0.91 0.65 11 0.0075 0.0054 0.62
## GAD_4 0.93 0.93 0.92 0.67 12 0.0068 0.0110 0.62
## GAD_5 0.93 0.93 0.93 0.70 14 0.0060 0.0081 0.71
## GAD_6 0.93 0.93 0.93 0.70 14 0.0060 0.0084 0.71
## GAD_7 0.93 0.93 0.93 0.68 13 0.0067 0.0100 0.66
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## GAD_1 275 0.90 0.89 0.88 0.86 1.01 0.99
## GAD_2 275 0.90 0.89 0.89 0.86 0.87 1.02
## GAD_3 275 0.91 0.90 0.91 0.87 0.99 1.05
## GAD_4 275 0.85 0.85 0.82 0.80 0.95 1.00
## GAD_5 275 0.77 0.78 0.73 0.70 0.51 0.88
## GAD_6 275 0.78 0.78 0.72 0.70 0.83 0.93
## GAD_7 275 0.84 0.84 0.80 0.78 0.71 0.93
##
## Non missing response frequency for each item
## 0 1 1.66666666666667 2 3 miss
## GAD_1 0.37 0.36 0 0.16 0.11 0
## GAD_2 0.47 0.30 0 0.11 0.12 0
## GAD_3 0.43 0.29 0 0.16 0.13 0
## GAD_4 0.41 0.33 0 0.14 0.11 0
## GAD_5 0.69 0.16 0 0.09 0.06 0
## GAD_6 0.45 0.33 0 0.13 0.08 0
## GAD_7 0.55 0.28 0 0.11 0.07 0
psych::alpha(qdf[, grep("RRS", names(qdf))])
##
## Reliability analysis
## Call: psych::alpha(x = qdf[, grep("RRS", names(qdf))])
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.9 0.9 0.92 0.46 8.6 0.0088 2 0.71 0.48
##
## 95% confidence boundaries
## lower alpha upper
## Feldt 0.88 0.9 0.92
## Duhachek 0.88 0.9 0.92
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## RRS_1 0.89 0.89 0.91 0.47 8.1 0.0093 0.029 0.50
## RRS_2 0.88 0.88 0.90 0.45 7.2 0.0101 0.029 0.44
## RRS_3 0.88 0.88 0.90 0.45 7.4 0.0101 0.028 0.44
## RRS_4 0.89 0.88 0.90 0.45 7.5 0.0098 0.028 0.48
## RRS_5 0.91 0.91 0.92 0.52 9.7 0.0085 0.014 0.51
## RRS_6 0.89 0.89 0.91 0.46 7.7 0.0097 0.028 0.47
## RRS_7 0.88 0.88 0.90 0.45 7.4 0.0101 0.026 0.47
## RRS_8 0.89 0.88 0.90 0.45 7.5 0.0100 0.026 0.46
## RRS_9 0.88 0.88 0.90 0.44 7.2 0.0103 0.029 0.44
## RRS_10 0.89 0.89 0.90 0.46 7.8 0.0094 0.030 0.49
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## RRS_1 275 0.67 0.66 0.62 0.58 1.8 0.94
## RRS_2 275 0.79 0.80 0.78 0.73 2.2 0.98
## RRS_3 275 0.78 0.77 0.75 0.72 2.0 0.98
## RRS_4 275 0.75 0.76 0.75 0.68 2.2 0.95
## RRS_5 275 0.39 0.43 0.33 0.30 1.4 0.70
## RRS_6 275 0.72 0.71 0.67 0.64 2.4 1.02
## RRS_7 275 0.79 0.77 0.76 0.72 2.0 1.07
## RRS_8 275 0.78 0.76 0.74 0.70 2.1 1.08
## RRS_9 275 0.81 0.81 0.79 0.75 2.0 1.03
## RRS_10 275 0.69 0.71 0.68 0.61 2.2 0.96
##
## Non missing response frequency for each item
## 1 1.33333333333333 2 2.22222222222222 3 4 miss
## RRS_1 0.51 0 0.30 0 0.11 0.08 0
## RRS_2 0.27 0 0.37 0 0.24 0.12 0
## RRS_3 0.40 0 0.32 0 0.19 0.09 0
## RRS_4 0.28 0 0.40 0 0.22 0.11 0
## RRS_5 0.67 0 0.27 0 0.04 0.03 0
## RRS_6 0.21 0 0.35 0 0.25 0.19 0
## RRS_7 0.44 0 0.27 0 0.15 0.14 0
## RRS_8 0.40 0 0.27 0 0.19 0.15 0
## RRS_9 0.40 0 0.29 0 0.18 0.12 0
## RRS_10 0.28 0 0.36 0 0.25 0.11 0
Sum scores
qdf$BDI_sum <- rowSums(qdf[, grep("BDI", names(qdf))])
qdf$GAD_sum <- rowSums(qdf[, grep("GAD", names(qdf))])
qdf$RRS_sum <- rowSums(qdf[, grep("RRS", names(qdf))])
qdf$depr_anx_sum <- qdf$BDI_sum + qdf$GAD_sum
#any(is.na(qdf))
#write.csv(qdf, "../data/questionnaire_df_with_impute.csv")
Put depr/anxisety sum and GAD sum into learn df
unique_ids <- unique(learn_df$ID)
for (i in 1:length(unique_ids)) {
learn_df[learn_df$ID==unique_ids[i], "depr_anx_sum"] <-
qdf[qdf$ID==unique_ids[i], "depr_anx_sum"]
learn_df[learn_df$ID==unique_ids[i], "bdi_sum"] <- qdf[qdf$ID==unique_ids[i], "BDI_sum"]
learn_df[learn_df$ID==unique_ids[i], "gad_sum"] <- qdf[qdf$ID==unique_ids[i], "GAD_sum"]
learn_df[learn_df$ID==unique_ids[i], "rrs_sum"] <- qdf[qdf$ID==unique_ids[i], "RRS_sum"]
test_df[test_df$ID==unique_ids[i], "depr_anx_sum"] <-
qdf[qdf$ID==unique_ids[i], "depr_anx_sum"]
test_df[test_df$ID==unique_ids[i], "bdi_sum"] <- qdf[qdf$ID==unique_ids[i], "BDI_sum"]
test_df[test_df$ID==unique_ids[i], "gad_sum"] <- qdf[qdf$ID==unique_ids[i], "GAD_sum"]
test_df[test_df$ID==unique_ids[i], "rrs_sum"] <- qdf[qdf$ID==unique_ids[i], "RRS_sum"]
}
# Check
# hist(learn_df$depr_anx_sum)
# hist(test_df$depr_anx_sum)
#
# hist(learn_df$rrs_sum)
# hist(test_df$rrs_sum)
# Export — df imported into for brms_s
# write.csv(learn_df, "../data/learn_df_with_qairre_data.csv")
# write.csv(test_df, "../data/test_df_with_qairre_data.csv")
Plot for supp fig
bdi_means <- qdf %>% group_by(group) %>% summarize(m=mean(BDI_sum))
bdi_p <- ggplot(qdf, aes(x=BDI_sum, fill=group, color=group)) +
geom_vline(data=bdi_means, aes(xintercept = m, color=group), size=3) +
geom_density(alpha=.4, linewidth=3) + ga + ap + lp + ylab("") +
xlab("") +
theme(axis.text.y = element_blank()) +
theme(axis.ticks.y = element_blank()) +
scale_fill_manual(values=c("red", "blue", "gray")) +
scale_color_manual(values=c("red", "blue", "gray")) +
ggtitle("Depression (BDI-II)") + tp +
tol
Just for legend
# ggplot(qdf, aes(x=BDI_sum, fill=group, color=group)) +
# geom_density(alpha=.4, linewidth=3) + ga + ap + lp + ylab("") +
# xlab("") +
# theme(axis.text.y = element_blank()) +
# theme(axis.ticks.y = element_blank()) +
# scale_fill_manual(values=c("red", "blue", "gray")) +
# scale_color_manual(values=c("red", "blue", "gray")) +
# ggtitle("Depression (BDI-II)") + tp
gad_means <- qdf %>% group_by(group) %>% summarize(m=mean(GAD_sum))
gad_p <- ggplot(qdf, aes(x=GAD_sum, fill=group, color=group)) +
geom_vline(data=gad_means, aes(xintercept = m, color=group), size=3) +
geom_density(alpha=.4, linewidth=3) + ga + ap + lp + ylab("") +
xlab("") +
theme(axis.text.y = element_blank()) +
theme(axis.ticks.y = element_blank()) +
scale_fill_manual(values=c("red", "blue", "gray")) +
scale_color_manual(values=c("red", "blue", "gray")) +
ggtitle("Generalized anxiety (GAD-7)") + tp +
tol
rrs_means <- qdf %>% group_by(group) %>% summarize(m=mean(RRS_sum))
rrs_p <- ggplot(qdf, aes(x=RRS_sum, fill=group, color=group)) +
geom_vline(data=rrs_means, aes(xintercept = m, color=group), size=3) +
geom_density(alpha=.4, linewidth=3) + ga + ap + lp + ylab("") +
xlab("") +
theme(axis.text.y = element_blank()) +
theme(axis.ticks.y = element_blank()) +
scale_fill_manual(values=c("red", "blue", "gray")) +
scale_color_manual(values=c("red", "blue", "gray")) +
ggtitle("Rumination (RRS-SF)") + tp +
tol
sum(table(qdf$group))
## [1] 275
sym_plot <- bdi_p/ gad_p / rrs_p
sym_plot
#ggsave("../paper/figs/supp-figs/symptom_plot.png", sym_plot, width=8.5, height = 8.5, dpi=200)
cor.test(learn_m$m, qdf$depr_anx_sum)
##
## Pearson's product-moment correlation
##
## data: learn_m$m and qdf$depr_anx_sum
## t = -0.01433, df = 273, p-value = 0.9886
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.1191390 0.1174287
## sample estimates:
## cor
## -0.0008672926
cor.test(test_m$m, qdf$depr_anx_sum)
##
## Pearson's product-moment correlation
##
## data: test_m$m and qdf$depr_anx_sum
## t = 0.14914, df = 273, p-value = 0.8816
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.1093747 0.1271742
## sample estimates:
## cor
## 0.009026022
cor.test(learn_m$m, qdf$RRS_sum)
##
## Pearson's product-moment correlation
##
## data: learn_m$m and qdf$RRS_sum
## t = -0.023301, df = 273, p-value = 0.9814
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.1196743 0.1168932
## sample estimates:
## cor
## -0.001410259
cor.test(test_m$m, qdf$RRS_sum)
##
## Pearson's product-moment correlation
##
## data: test_m$m and qdf$RRS_sum
## t = 0.0072235, df = 273, p-value = 0.9942
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.1178529 0.1187150
## sample estimates:
## cor
## 0.0004371846
ctest <- cor.test(qdf$GAD_sum, qdf$BDI_sum)
r <- round(ctest$estimate, 2)
p <- round(ctest$p.value, 2)
str <- paste("r =", r, "p =", p)
ggplot(qdf, aes(x=GAD_sum, BDI_sum)) +
geom_smooth(method='lm', formula= y~x, color="black", size=3) +
geom_point(size=6, fill="white", pch=21) +
ga + ap + tp +
ggtitle("", subtitle=str) +
ylab("BDI-II") + xlab("GAD-7") +
theme(plot.subtitle = element_text(size = 25, face = "bold"))
Create high low DA, BDI and GAD
median(qdf$BDI_sum)
## [1] 10
median(qdf$GAD_sum)
## [1] 4
hist(qdf$depr_anx_sum, breaks=100)
median(qdf$depr_anx_sum)
## [1] 15
table(if_else(qdf$depr_anx_sum > median(qdf$depr_anx_sum), 1, 0))
##
## 0 1
## 138 137
#table(if_else(qdf$depr_anx_sum > median(qdf$depr_anx_sum), 1, 0)) # Check that's about half
qdf$DA_hilo <- if_else(qdf$depr_anx_sum > median(qdf$depr_anx_sum), 1, 0)
qdf %>% group_by(DA_hilo) %>% summarize(m=median(depr_anx_sum)) # Confirm
## # A tibble: 2 × 2
## DA_hilo m
## <dbl> <dbl>
## 1 0 4
## 2 1 33
#table(if_else(qdf$depr_anx_sum > median(qdf$depr_anx_sum), 1, 0))
qdf$dep_hilo <- if_else(qdf$BDI_sum > median(qdf$BDI_sum), 1, 0)
#qdf %>% group_by(dep_hilo) %>% summarize(m=median(depr_anx_sum))
#table(if_else(qdf$GAD_sum > median(qdf$GAD_sum), 1, 0))
qdf$anx_hilo <- if_else(qdf$GAD_sum > median(qdf$GAD_sum), 1, 0)
#qdf %>% group_by(anx_hilo) %>% summarize(m=median(GAD_sum))
table(if_else(qdf$RRS_sum > median(qdf$RRS_sum), 1, 0))
##
## 0 1
## 146 129
qdf$rrs_hilo <- if_else(qdf$RRS_sum > median(qdf$RRS_sum), 1, 0)
qdf %>% group_by(rrs_hilo) %>% summarize(m=median(RRS_sum))
## # A tibble: 2 × 2
## rrs_hilo m
## <dbl> <dbl>
## 1 0 15
## 2 1 25
Put into learning and test dfs
high_depr_anx_pts <- qdf[qdf$DA_hilo == 1, "ID"]
low_depr_anx_pts <- qdf[qdf$DA_hilo == 0, "ID"]
learn_df[learn_df$ID %in% high_depr_anx_pts, "DA_hilo"] <- "high depr-anx"
learn_df[learn_df$ID %in% low_depr_anx_pts, "DA_hilo"] <- "low depr-anx"
test_df[test_df$ID %in% high_depr_anx_pts, "DA_hilo"] <- "high depr-anx"
test_df[test_df$ID %in% low_depr_anx_pts, "DA_hilo"] <- "low depr-anx"
high_rum_pts <- qdf[qdf$rrs_hilo == 1, "ID"]
low_rum_pts <- qdf[qdf$rrs_hilo == 0, "ID"]
learn_df[learn_df$ID %in% high_rum_pts, "rrs_hilo"] <- "high rum"
learn_df[learn_df$ID %in% low_rum_pts, "rrs_hilo"] <- "low rum"
test_df[test_df$ID %in% high_rum_pts, "rrs_hilo"] <- "high rum"
test_df[test_df$ID %in% low_rum_pts, "rrs_hilo"] <- "low rum"
Plot proportion correct by median split
pcor_ss_DA <- data.frame(learn_df %>% group_by(stim_iter, set_size,
ID, DA_hilo) %>% summarize(m=mean(correct)))
## `summarise()` has grouped output by 'stim_iter', 'set_size', 'ID'. You can
## override using the `.groups` argument.
pcor_ss_err_DA <- Rmisc::summarySEwithin(pcor_ss_DA,
measurevar = "m",
withinvars = c("set_size", "stim_iter", "DA_hilo"),
idvar = "ID")
## Automatically converting the following non-factors to factors: set_size, stim_iter, DA_hilo
pcor_ss_err_DA$DA_hilo <- factor(pcor_ss_err_DA$DA_hilo, levels=c("low depr-anx", "high depr-anx"))
pcor_by_grp <- ggplot(pcor_ss_err_DA, aes(x=stim_iter, y=m, group=as.factor(DA_hilo),
color=as.factor(DA_hilo))) +
geom_line() +
geom_ribbon(aes(ymin=m-se, ymax=m+se), fill='gray57', alpha=.45) +
geom_hline(yintercept=.33, size=1.5, color="gray57") + # chance line
geom_hline(yintercept=c(.5, .6, .7, .8, .9, 1), linetype="dotted") +
geom_point(aes(fill=as.factor(DA_hilo)), color="black", size=5, pch=21) +
geom_vline(xintercept=c(2, 5, 8, 10), linetype="dotted") +
#geom_point(aes(fill=as.factor(set_size)), color="black", size=5, pch=21) +
annotate("rect", xmin=6, xmax=10.5, ymin=.3, ymax=1.1, alpha=0.2, fill="gray57") +
ga + ap + lp + xlab("Stimulus iteration") + ylab("Proportion correct") +
ggtitle("") + tp + facet_wrap(~ set_size) + ft + lp +
scale_fill_manual(values=c("orange1", "brown")) +
scale_color_manual(values=c("orange1", "brown")) + tol
pcor_by_grp
#ggsave("../paper/figs/pieces/pcor_by_grp.png", pcor_by_grp, width=14, height=10, dpi=300)
pcor_ss_rrs <- data.frame(learn_df %>% group_by(stim_iter, set_size,
ID, rrs_hilo) %>% summarize(m=mean(correct)))
## `summarise()` has grouped output by 'stim_iter', 'set_size', 'ID'. You can
## override using the `.groups` argument.
pcor_ss_err_rrs <- Rmisc::summarySEwithin(pcor_ss_rrs,
measurevar = "m",
withinvars = c("set_size", "stim_iter", "rrs_hilo"),
idvar = "ID")
## Automatically converting the following non-factors to factors: set_size, stim_iter, rrs_hilo
pcor_ss_err_rrs$rrs_hilo <- factor(pcor_ss_err_rrs$rrs_hilo, levels=c("low rum", "high rum"))
pcor_by_grp_rum <- ggplot(pcor_ss_err_rrs, aes(x=stim_iter, y=m, group=as.factor(rrs_hilo),
color=as.factor(rrs_hilo))) +
geom_line() +
geom_ribbon(aes(ymin=m-se, ymax=m+se), fill='gray57', alpha=.45) +
geom_hline(yintercept=.33, size=1.5, color="gray57") + # chance line
geom_hline(yintercept=c(.5, .6, .7, .8, .9, 1), linetype="dotted") +
geom_point(aes(fill=as.factor(rrs_hilo)), color="black", size=5, pch=21) +
geom_vline(xintercept=c(2, 5, 8, 10), linetype="dotted") +
#geom_point(aes(fill=as.factor(set_size)), color="black", size=5, pch=21) +
annotate("rect", xmin=6, xmax=10.5, ymin=.3, ymax=1.1, alpha=0.2, fill="gray57") +
ga + ap + lp + xlab("Stimulus iteration") + ylab("Proportion correct") +
ggtitle("") + tp + facet_wrap(~ set_size) + ft + lp +
scale_fill_manual(values=c("red", "darkred")) +
scale_color_manual(values=c("red", "darkred"))
Rumination version
pcor_by_grp_rum
#ggsave("../paper/figs/pieces/pcor_by_grp_rum.png", pcor_by_grp_rum, width=14, height=8, dpi=300)
test_error_df <- test_df %>% filter(correct==0)
np_test_summs_da <- test_error_df %>% group_by(ID, DA_hilo, phase) %>%
summarize(mw=mean(worst), mn=mean(neutral), n=n()) %>% mutate(neutral_pref=(mn-mw))
## `summarise()` has grouped output by 'ID', 'DA_hilo'. You can override using the
## `.groups` argument.
np_test_si_da <- Rmisc::summarySEwithin(np_test_summs_da,
measurevar = "neutral_pref",
withinvars = c("DA_hilo", "phase"),
idvar = "ID")
## Automatically converting the following non-factors to factors: DA_hilo, phase
np_test_si_da$DA_hilo <- factor(np_test_si_da$DA_hilo, levels=c("low depr-anx", "high depr-anx"))
ggplot(np_test_si_da, aes(x=phase, y=neutral_pref, fill=phase)) +
geom_hline(yintercept=.0, size=1.5, color="gray57") + # chance line
geom_bar(stat="identity", color="black") +
geom_errorbar(aes(ymin=neutral_pref-se, ymax=neutral_pref+se), width=.2) +
ga + ap + tol + xlab("Phase") + ylab("") + tp +
ggtitle(#"Neutral preference during test
"Empirical") + scale_fill_manual(values=c("gray81", "gray40")) + ylim(-.05, .23) + facet_wrap(~ DA_hilo) + ft
brms_s.R for
behaviorPerformace in terms of proprortion correct
phase1_learn_perf <- read.csv("../model_res/brms_res/reduced_perf_model_phase1__86050__.csv")
phase2_learn_perf <- read.csv("../model_res/brms_res/reduced_correct_perf_model_phase2__76459__.csv")
test_perf <- read.csv("../model_res/brms_res/reduced_correct_test_model__24237__.csv")
ReturnPosteriorMeanAnd90CI <- function(posterior) {
cat("\nposterior mean =", mean(posterior))
cat("\n90% CI ="); print (bayestestR::ci(posterior, ci = .9, method = "HDI"))
}
Neutral preference
np_da <- read.csv("../model_res/brms_res/neutral_pref_learn_model__64424__.csv")
# Version that converged that has intercept set to 0 but models REs in depr/anx
np_test_da <- read.csv("../model_res/brms_res/neutral_pref_test_model_0int__14655__.csv")
Rumination
phase1_rum <- read.csv("../model_res/brms_res/RUM_reduced_perf_model_phase1__63279__.csv")
phase2_rum <- read.csv("../model_res/brms_res/RUM_reduced_correct_perf_model_phase2__25543__.csv")
test_rum <- read.csv("../model_res/brms_res/RUM_reduced_correct_test_model__53350__.csv")
np_rum <- read.csv("../model_res/brms_res/RUM_neutral_pref_learn_model__89767__.csv")
np_test_rum <- read.csv("../model_res/brms_res/RUM_neutral_pref_test_model__71650__.csv")
Phase 1
cat("\n Depression prop correct phase 1: \n")
##
## Depression prop correct phase 1:
# DA
ReturnPosteriorMeanAnd90CI(phase1_learn_perf$b_scaledepr_anx_sum)
##
## posterior mean = -0.005247766
## 90% CI =90% HDI: [-0.04, 0.03]
ReturnPosteriorMeanAnd90CI(phase1_learn_perf$b_scaleset_size) # Sanity check this is negative and tight
##
## posterior mean = -0.3166876
## 90% CI =90% HDI: [-0.33, -0.30]
ReturnPosteriorMeanAnd90CI(phase1_learn_perf$b_scaledepr_anx_sum.scaleset_size)
##
## posterior mean = 0.004837842
## 90% CI =90% HDI: [-0.01, 0.02]
cat("\n\n Rumination prop correct phase 1: \n")
##
##
## Rumination prop correct phase 1:
# Rum
ReturnPosteriorMeanAnd90CI(phase1_rum$b_scalerrs_sum)
##
## posterior mean = -0.007088953
## 90% CI =90% HDI: [-0.04, 0.03]
ReturnPosteriorMeanAnd90CI(phase1_rum$b_scaleset_size) # Sanity check - v consistent with above
##
## posterior mean = -0.3169965
## 90% CI =90% HDI: [-0.33, -0.30]
ReturnPosteriorMeanAnd90CI(phase1_rum$b_scalerrs_sum.scaleset_size)
##
## posterior mean = 0.01584776
## 90% CI =90% HDI: [ 0.00, 0.03]
Phase 2
# DA
cat("\n\n Depression prop correct phase 2: \n")
##
##
## Depression prop correct phase 2:
ReturnPosteriorMeanAnd90CI(phase2_learn_perf$b_scaledepr_anx_sum)
##
## posterior mean = 6.641381e-05
## 90% CI =90% HDI: [-0.08, 0.07]
ReturnPosteriorMeanAnd90CI(phase2_learn_perf$b_scaleset_size) # Sanity check this is negative and tight
##
## posterior mean = -0.2150235
## 90% CI =90% HDI: [-0.24, -0.19]
ReturnPosteriorMeanAnd90CI(phase2_learn_perf$b_scaledepr_anx_sum.scaleset_size)
##
## posterior mean = 0.02980973
## 90% CI =90% HDI: [0.01, 0.05]
# Rum
cat("\n\n Rumination prop correct phase 2: \n")
##
##
## Rumination prop correct phase 2:
ReturnPosteriorMeanAnd90CI(phase2_rum$b_scalerrs_sum)
##
## posterior mean = 0.0003879612
## 90% CI =90% HDI: [-0.08, 0.08]
ReturnPosteriorMeanAnd90CI(phase2_rum$b_scaleset_size) # Sanity check - v consistent with above
##
## posterior mean = -0.2149793
## 90% CI =90% HDI: [-0.24, -0.19]
ReturnPosteriorMeanAnd90CI(phase2_rum$b_scalerrs_sum.scaleset_size)
##
## posterior mean = 0.01623668
## 90% CI =90% HDI: [-0.01, 0.04]
Test
# DA
cat("\n\n Depression prop correct test: \n")
##
##
## Depression prop correct test:
ReturnPosteriorMeanAnd90CI(test_perf$b_scaledepr_anx_sum)
##
## posterior mean = 0.009176028
## 90% CI =90% HDI: [-0.08, 0.09]
# Rum
cat("\n\n Rumination prop correct test: \n")
##
##
## Rumination prop correct test:
ReturnPosteriorMeanAnd90CI(test_rum$b_scalerrs_sum)
##
## posterior mean = 0.004301401
## 90% CI =90% HDI: [-0.08, 0.09]
The intercept posteriors show that the evidence for a neutral preference found in the frequentist models holds in these Bayesian regression models, now while also controlling for psychiatric symptoms
cat("\n Depression : \n")
##
## Depression :
cat("\n Evidence controlling for depression : "); ReturnPosteriorMeanAnd90CI(np_da$b_Intercept)
##
## Evidence controlling for depression :
##
## posterior mean = 0.08774017
## 90% CI =90% HDI: [0.06, 0.11]
cat("\n Evidence for depession effect : ");ReturnPosteriorMeanAnd90CI(np_da$b_scaledepr_anx_sum)
##
## Evidence for depession effect :
##
## posterior mean = -0.01066593
## 90% CI =90% HDI: [-0.03, 0.01]
cat("\n Rumination : \n")
##
## Rumination :
cat("\n Evidence controlling for rumination : "); ReturnPosteriorMeanAnd90CI(np_rum$b_Intercept)
##
## Evidence controlling for rumination :
##
## posterior mean = 0.08747077
## 90% CI =90% HDI: [0.06, 0.11]
cat("\n Evidence for rumination effect : ");ReturnPosteriorMeanAnd90CI(np_rum$b_scalerrs_sum)
##
## Evidence for rumination effect :
##
## posterior mean = 0.004026421
## 90% CI =90% HDI: [-0.02, 0.03]
Note that HDI for intercept (indicating neutral pref at test) includes 0 (although it’s close controlling for depression symptoms)
cat("\n Depression : \n")
##
## Depression :
cat("\n Evidence controlling for depression : "); ReturnPosteriorMeanAnd90CI(np_test_da$b_Intercept)
##
## Evidence controlling for depression :
##
## posterior mean = 0.03915869
## 90% CI =90% HDI: [ 0.00, 0.09]
cat("\n Evidence for depession effect : ");ReturnPosteriorMeanAnd90CI(np_test_da$b_scaledepr_anx_sum)
##
## Evidence for depession effect :
##
## posterior mean = -0.007424603
## 90% CI =90% HDI: [-0.06, 0.05]
cat("\n Rumination : \n")
##
## Rumination :
cat("\n Evidence controlling for rumination : "); ReturnPosteriorMeanAnd90CI(np_test_rum$b_Intercept)
##
## Evidence controlling for rumination :
##
## posterior mean = 0.03967129
## 90% CI =90% HDI: [-0.01, 0.09]
cat("\n Evidence for rumination effect : ");ReturnPosteriorMeanAnd90CI(np_test_rum$b_scalerrs_sum)
##
## Evidence for rumination effect :
##
## posterior mean = -0.00610522
## 90% CI =90% HDI: [-0.06, 0.05]
Make sure Rhats below 1.1
CheckRhat("../model_res/brms_res/Rhatreduced_correct_perf_model_phase2__76459__.csv")
## [1] TRUE
CheckRhat("../model_res/brms_res/Rhatreduced_perf_model_phase1__86050__.csv")
## [1] TRUE
CheckRhat("../model_res/brms_res/Rhatreduced_correct_test_model__24237__.csv")
## [1] TRUE
CheckRhat("../model_res/brms_res/Rhatneutral_pref_learn_model__64424__.csv")
## [1] TRUE
CheckRhat("../model_res/brms_res/Rhatneutral_pref_test_model__70753__.csv")
## [1] TRUE
CheckRhat("../model_res/brms_res/RhatRUM_neutral_pref_learn_model__89767__.csv")
## [1] TRUE
CheckRhat("../model_res/brms_res/RhatRUM_neutral_pref_test_model__71650__.csv")
## [1] TRUE
CheckRhat("../model_res/brms_res/RhatRUM_reduced_correct_perf_model_phase2__25543__.csv")
## [1] TRUE
CheckRhat("../model_res/brms_res/RhatRUM_reduced_perf_model_phase1__63279__.csv")
## [1] TRUE
CheckRhat("../model_res/brms_res/RhatRUM_reduced_correct_test_model__53350__.csv")
## [1] TRUE
Get ready to plot proportion correct posterior for depr/anxiety …
depr_anx_learn_perf_posteriors <- rbind(
data.frame("trace"=phase1_learn_perf$b_scaledepr_anx_sum, "label"="da_main_effect", "phase"="Learn Phase 1"),
data.frame("trace"=phase2_learn_perf$b_scaledepr_anx_sum, "label"="da_main_effect", "phase"="Learn Phase 2"),
data.frame("trace"=phase1_learn_perf$b_scaledepr_anx_sum.scaleset_size,
"label"="da_ss_int", "phase"="Learn Phase 1"),
data.frame("trace"=phase2_learn_perf$b_scaledepr_anx_sum.scaleset_size,
"label"="da_ss_int", "phase"="Learn Phase 2")
)
depr_anx_test_perf_posterior <- rbind(
data.frame("trace"=test_perf$b_scaledepr_anx_sum, "label"="da_main_effect", "phase"="Test")
)
learn_correct_posterior_p <- ggplot(depr_anx_learn_perf_posteriors,
aes(x=trace, color=label)) +
geom_vline(xintercept=0, color="gray20", size=1.5) +
geom_density(fill="gray57", alpha=.15, size=1.5) +
facet_wrap(~ phase, ncol=1) + ga + ap +
scale_color_manual(values=c("tan1", "chocolate")) + ylab("") +
theme(axis.text.x = element_text(size=15), axis.ticks=element_blank(), axis.text.y=element_blank()) +
#labs(title="Learning") +
labs(title="") +
tol + xlab("") + theme(strip.text.x = element_text(size = 15)) + tp
learn_correct_posterior_p
test_correct_posterior_p <-ggplot(depr_anx_test_perf_posterior,
aes(x=trace, color=label)) +
geom_vline(xintercept=0, color="gray20", size=1.5) +
geom_density(fill="gray57", alpha=.15, size=1.5) +
facet_wrap(~ phase, ncol=1) +
ga + ap +
scale_color_manual(values=c("tan1")) +
ylab("") +
theme(axis.text.x = element_text(size=15),
axis.ticks=element_blank(), axis.text.y=element_blank()) +
#labs(title="Test") +
tol + xlab("") + theme(strip.text.x = element_text(size = 15)) + tp #+
test_correct_posterior_p
# ggsave("../paper/figs/pieces/fig7_learn_correct_posterior_p.png", learn_correct_posterior_p, width=2.2, height=3.3, dpi=500)
# ggsave("../paper/figs/pieces/fig7_test_correct_posterior_p.png",
# test_correct_posterior_p, width=2.2, height=1.65, dpi=500)
… and rum
rum_learn_perf_posteriors <- rbind(
data.frame("trace"=phase1_rum$b_scalerrs_sum, "label"="rum_main_effect", "phase"="Learn Phase 1"),
data.frame("trace"=phase2_rum$b_scalerrs_sum, "label"="rum_main_effect", "phase"="Learn Phase 2"),
data.frame("trace"=phase1_rum$b_scalerrs_sum.scaleset_size,
"label"="rum_ss_int", "phase"="Learn Phase 1"),
data.frame("trace"=phase2_rum$b_scalerrs_sum.scaleset_size,
"label"="rum_ss_int", "phase"="Learn Phase 2")
)
rum_test_perf_posterior <- rbind(
data.frame("trace"=test_rum$b_scalerrs_sum, "label"="rum_main_effect", "phase"="Test")
)
learn_correct_posterior_p_RUM <-
ggplot(rum_learn_perf_posteriors,
aes(x=trace, color=label)) +
geom_vline(xintercept=0, color="gray20", size=1.5) +
geom_density(fill="gray57", alpha=.15, size=1.5) +
facet_wrap(~ phase, ncol=1) + ga + ap +
scale_color_manual(values=c("indianred2", "darkred")) + ylab("") +
theme(axis.text.x = element_text(size=15), axis.ticks=element_blank(), axis.text.y=element_blank()) +
#labs(title="Learning") +
labs(title="") +
tol + xlab("") + theme(strip.text.x = element_text(size = 15)) + tp
learn_correct_posterior_p_RUM
test_correct_posterior_p_RUM <-
ggplot(rum_test_perf_posterior,
aes(x=trace, color=label)) +
geom_vline(xintercept=0, color="gray20", size=1.5) +
geom_density(fill="gray57", alpha=.15, size=1.5) +
facet_wrap(~ phase, ncol=1) +
ga + ap +
scale_color_manual(values=c("indianred1")) +
ylab("") +
theme(axis.text.x = element_text(size=15), axis.ticks=element_blank(), axis.text.y=element_blank()) +
#labs(title="Test") +
tol + xlab("") + theme(strip.text.x = element_text(size = 15)) + tp #+
test_correct_posterior_p_RUM
# ggsave("../paper/figs/pieces/fig7_learn_correct_posterior_p_RUM.png", learn_correct_posterior_p_RUM, width=2.2, height=3.3, dpi=500)
# ggsave("../paper/figs/pieces/fig7_test_correct_posterior_p_RUM.png",
# test_correct_posterior_p_RUM, width=2.2, height=1.65, dpi=500)
Get ready to plot neutral pref posterior for depr/anxiety and rum
depr_anx_np_posteriors <- rbind(
data.frame("trace"=np_da$b_scaledepr_anx_sum, "label"="da_effect", "phase"="Learn"),
data.frame("trace"=np_test_da$b_scaledepr_anx_sum, "label"="da_main_effect", "phase"="Test")
)
rum_np_posteriors <- rbind(
data.frame("trace"=np_rum$b_scalerrs_sum, "label"="rum_effect", "phase"="Learn"),
data.frame("trace"=np_test_rum$b_scalerrs_sum, "label"="da_main_effect", "phase"="Test")
)
np_da_posterior_p <-
ggplot(depr_anx_np_posteriors,
aes(x=trace, color=label)) +
geom_vline(xintercept=0, color="gray20", size=1.5) +
geom_density(fill="gray57", alpha=.15, size=1.5) +
facet_wrap(~ phase, ncol=1) + ga + ap +
scale_color_manual(values=c("tan1", "tan1")) + ylab("") +
theme(axis.text.x = element_text(size=10), axis.ticks=element_blank(), axis.text.y=element_blank()) +
#labs(title="Learning") +
labs(title="") +
tol + xlab("") + theme(strip.text.x = element_text(size = 15)) + tp
np_da_posterior_p
rum_np_posterior_p <-
ggplot(rum_np_posteriors,
aes(x=trace, color=label)) +
geom_vline(xintercept=0, color="gray20", size=1.5) +
geom_density(fill="gray57", alpha=.15, size=1.5) +
facet_wrap(~ phase, ncol=1) + ga + ap +
scale_color_manual(values=c("indianred2", "indianred2")) + ylab("") +
theme(axis.text.x = element_text(size=10), axis.ticks=element_blank(), axis.text.y=element_blank()) +
#labs(title="Learning") +
labs(title="") +
tol + xlab("") + theme(strip.text.x = element_text(size = 15)) + tp
rum_np_posterior_p
# ggsave("../paper/figs/pieces/fig7_np_da_posterior_p.png", np_da_posterior_p ,
# width=2.2, height=3.3, dpi=500)
# ggsave("../paper/figs/pieces/fig7_rum_np_posterior_p.png",
# rum_np_posterior_p, width=2.2, height=3.3, dpi=500)
learn_m <- learn_df %>% group_by(ID) %>% summarize(m=mean(correct))
assert(all(learn_m$ID==qdf$ID))
test_m <- test_df %>% group_by(ID) %>% summarize(m=mean(correct))
assert(all(learn_m$ID==qdf$ID))
assert(all(m35$ID==qdf$ID))
assert(all(m35$ID==learn_m$ID))
assert(all(m35$ID==test_m$ID))
m35$depr_anx_sum <- qdf$depr_anx_sum
m35$BDI_sum <- qdf$BDI_sum
m35$GAD_sum <- qdf$GAD_sum
m35$RRS_sum <- qdf$RRS_sum
m35$learn_m <- learn_m$m
m35$test_m <- test_m$m
Many parameters predict learning and test performance…
summary(learn_preds_bayes <- rstanarm::stan_glm(
scale(learn_m) ~
scale(kappa) +
scale(alpha_pos) +
scale(alpha_neg) +
scale(phi) +
scale(rho) +
scale(rl_off) +
scale(epsilon) +
scale(not_pun_bonus),
data=m35))
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 4.1e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.41 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.034 seconds (Warm-up)
## Chain 1: 0.052 seconds (Sampling)
## Chain 1: 0.086 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 6e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.032 seconds (Warm-up)
## Chain 2: 0.054 seconds (Sampling)
## Chain 2: 0.086 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 9e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.09 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.031 seconds (Warm-up)
## Chain 3: 0.052 seconds (Sampling)
## Chain 3: 0.083 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 9e-06 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.09 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.033 seconds (Warm-up)
## Chain 4: 0.047 seconds (Sampling)
## Chain 4: 0.08 seconds (Total)
## Chain 4:
##
## Model Info:
## function: stan_glm
## family: gaussian [identity]
## formula: scale(learn_m) ~ scale(kappa) + scale(alpha_pos) + scale(alpha_neg) +
## scale(phi) + scale(rho) + scale(rl_off) + scale(epsilon) +
## scale(not_pun_bonus)
## algorithm: sampling
## sample: 4000 (posterior sample size)
## priors: see help('prior_summary')
## observations: 275
## predictors: 9
##
## Estimates:
## mean sd 10% 50% 90%
## (Intercept) 0.0 0.0 0.0 0.0 0.0
## scale(kappa) 0.3 0.0 0.2 0.3 0.3
## scale(alpha_pos) 0.4 0.0 0.4 0.4 0.5
## scale(alpha_neg) 0.2 0.0 0.1 0.2 0.2
## scale(phi) -0.3 0.0 -0.3 -0.3 -0.2
## scale(rho) 0.1 0.0 0.1 0.1 0.2
## scale(rl_off) -0.3 0.0 -0.3 -0.3 -0.2
## scale(epsilon) -0.3 0.0 -0.4 -0.3 -0.3
## scale(not_pun_bonus) -0.1 0.0 -0.2 -0.1 -0.1
## sigma 0.5 0.0 0.5 0.5 0.6
##
## Fit Diagnostics:
## mean sd 10% 50% 90%
## mean_PPD 0.0 0.0 -0.1 0.0 0.1
##
## The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
##
## MCMC diagnostics
## mcse Rhat n_eff
## (Intercept) 0.0 1.0 4451
## scale(kappa) 0.0 1.0 3653
## scale(alpha_pos) 0.0 1.0 3728
## scale(alpha_neg) 0.0 1.0 4008
## scale(phi) 0.0 1.0 3418
## scale(rho) 0.0 1.0 3689
## scale(rl_off) 0.0 1.0 4010
## scale(epsilon) 0.0 1.0 3723
## scale(not_pun_bonus) 0.0 1.0 3857
## sigma 0.0 1.0 3701
## mean_PPD 0.0 1.0 4174
## log-posterior 0.1 1.0 1740
##
## For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).
summary(test_preds_bayes <- rstanarm::stan_glm(
scale(test_m) ~
scale(kappa) +
scale(alpha_pos) +
scale(alpha_neg) +
scale(phi) +
scale(rho) +
scale(rl_off) +
scale(epsilon) +
scale(not_pun_bonus),
data=m35))
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 1.4e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.14 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.032 seconds (Warm-up)
## Chain 1: 0.048 seconds (Sampling)
## Chain 1: 0.08 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 7e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.07 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.033 seconds (Warm-up)
## Chain 2: 0.047 seconds (Sampling)
## Chain 2: 0.08 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 6e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.033 seconds (Warm-up)
## Chain 3: 0.044 seconds (Sampling)
## Chain 3: 0.077 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 6e-06 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.034 seconds (Warm-up)
## Chain 4: 0.053 seconds (Sampling)
## Chain 4: 0.087 seconds (Total)
## Chain 4:
##
## Model Info:
## function: stan_glm
## family: gaussian [identity]
## formula: scale(test_m) ~ scale(kappa) + scale(alpha_pos) + scale(alpha_neg) +
## scale(phi) + scale(rho) + scale(rl_off) + scale(epsilon) +
## scale(not_pun_bonus)
## algorithm: sampling
## sample: 4000 (posterior sample size)
## priors: see help('prior_summary')
## observations: 275
## predictors: 9
##
## Estimates:
## mean sd 10% 50% 90%
## (Intercept) 0.0 0.0 0.0 0.0 0.0
## scale(kappa) 0.1 0.0 0.0 0.1 0.1
## scale(alpha_pos) 0.6 0.0 0.6 0.6 0.7
## scale(alpha_neg) 0.2 0.0 0.1 0.2 0.2
## scale(phi) -0.2 0.0 -0.2 -0.2 -0.1
## scale(rho) 0.1 0.0 0.0 0.1 0.1
## scale(rl_off) -0.5 0.0 -0.5 -0.5 -0.4
## scale(epsilon) -0.2 0.0 -0.2 -0.2 -0.1
## scale(not_pun_bonus) 0.0 0.0 -0.1 0.0 0.0
## sigma 0.5 0.0 0.5 0.5 0.6
##
## Fit Diagnostics:
## mean sd 10% 50% 90%
## mean_PPD 0.0 0.0 -0.1 0.0 0.1
##
## The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
##
## MCMC diagnostics
## mcse Rhat n_eff
## (Intercept) 0.0 1.0 3958
## scale(kappa) 0.0 1.0 3307
## scale(alpha_pos) 0.0 1.0 3426
## scale(alpha_neg) 0.0 1.0 4644
## scale(phi) 0.0 1.0 3455
## scale(rho) 0.0 1.0 3955
## scale(rl_off) 0.0 1.0 3959
## scale(epsilon) 0.0 1.0 4398
## scale(not_pun_bonus) 0.0 1.0 4318
## sigma 0.0 1.0 4648
## mean_PPD 0.0 1.0 4140
## log-posterior 0.1 1.0 1777
##
## For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).
summary(da_preds_bayes <- rstanarm::stan_glm(
scale(depr_anx_sum) ~
scale(kappa) +
scale(alpha_pos) +
scale(alpha_neg) +
scale(phi) +
scale(rho) +
scale(rl_off) +
scale(epsilon) +
scale(not_pun_bonus),
data=m35))
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 1.4e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.14 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.03 seconds (Warm-up)
## Chain 1: 0.047 seconds (Sampling)
## Chain 1: 0.077 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 7e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.07 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.031 seconds (Warm-up)
## Chain 2: 0.049 seconds (Sampling)
## Chain 2: 0.08 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 7e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.07 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.059 seconds (Warm-up)
## Chain 3: 0.05 seconds (Sampling)
## Chain 3: 0.109 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 8e-06 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.08 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.032 seconds (Warm-up)
## Chain 4: 0.056 seconds (Sampling)
## Chain 4: 0.088 seconds (Total)
## Chain 4:
##
## Model Info:
## function: stan_glm
## family: gaussian [identity]
## formula: scale(depr_anx_sum) ~ scale(kappa) + scale(alpha_pos) + scale(alpha_neg) +
## scale(phi) + scale(rho) + scale(rl_off) + scale(epsilon) +
## scale(not_pun_bonus)
## algorithm: sampling
## sample: 4000 (posterior sample size)
## priors: see help('prior_summary')
## observations: 275
## predictors: 9
##
## Estimates:
## mean sd 10% 50% 90%
## (Intercept) 0.0 0.1 -0.1 0.0 0.1
## scale(kappa) -0.1 0.1 -0.2 -0.1 0.0
## scale(alpha_pos) 0.1 0.1 0.0 0.1 0.1
## scale(alpha_neg) 0.0 0.1 -0.1 0.0 0.1
## scale(phi) -0.1 0.1 -0.2 -0.1 0.0
## scale(rho) 0.0 0.1 -0.1 0.0 0.1
## scale(rl_off) 0.0 0.1 -0.1 0.0 0.1
## scale(epsilon) 0.0 0.1 0.0 0.0 0.1
## scale(not_pun_bonus) 0.1 0.1 0.0 0.1 0.1
## sigma 1.0 0.0 1.0 1.0 1.1
##
## Fit Diagnostics:
## mean sd 10% 50% 90%
## mean_PPD 0.0 0.1 -0.1 0.0 0.1
##
## The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
##
## MCMC diagnostics
## mcse Rhat n_eff
## (Intercept) 0.0 1.0 3809
## scale(kappa) 0.0 1.0 3137
## scale(alpha_pos) 0.0 1.0 3216
## scale(alpha_neg) 0.0 1.0 3906
## scale(phi) 0.0 1.0 3193
## scale(rho) 0.0 1.0 3812
## scale(rl_off) 0.0 1.0 4305
## scale(epsilon) 0.0 1.0 4081
## scale(not_pun_bonus) 0.0 1.0 4331
## sigma 0.0 1.0 3725
## mean_PPD 0.0 1.0 3910
## log-posterior 0.1 1.0 1725
##
## For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).
summary(rrs_preds_bayes <- rstanarm::stan_glm(
scale(RRS_sum) ~
scale(kappa) +
scale(alpha_pos) +
scale(alpha_neg) +
scale(phi) +
scale(rho) +
scale(rl_off) +
scale(epsilon) +
scale(not_pun_bonus),
data=m35))
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 1.1e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.11 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.031 seconds (Warm-up)
## Chain 1: 0.05 seconds (Sampling)
## Chain 1: 0.081 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 6e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.03 seconds (Warm-up)
## Chain 2: 0.047 seconds (Sampling)
## Chain 2: 0.077 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 6e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.032 seconds (Warm-up)
## Chain 3: 0.046 seconds (Sampling)
## Chain 3: 0.078 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 5e-06 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.05 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.034 seconds (Warm-up)
## Chain 4: 0.05 seconds (Sampling)
## Chain 4: 0.084 seconds (Total)
## Chain 4:
##
## Model Info:
## function: stan_glm
## family: gaussian [identity]
## formula: scale(RRS_sum) ~ scale(kappa) + scale(alpha_pos) + scale(alpha_neg) +
## scale(phi) + scale(rho) + scale(rl_off) + scale(epsilon) +
## scale(not_pun_bonus)
## algorithm: sampling
## sample: 4000 (posterior sample size)
## priors: see help('prior_summary')
## observations: 275
## predictors: 9
##
## Estimates:
## mean sd 10% 50% 90%
## (Intercept) 0.0 0.1 -0.1 0.0 0.1
## scale(kappa) 0.0 0.1 -0.1 0.0 0.1
## scale(alpha_pos) 0.0 0.1 -0.1 0.0 0.1
## scale(alpha_neg) 0.0 0.1 -0.1 0.0 0.1
## scale(phi) 0.0 0.1 -0.1 0.0 0.1
## scale(rho) 0.0 0.1 0.0 0.0 0.1
## scale(rl_off) 0.0 0.1 0.0 0.0 0.1
## scale(epsilon) 0.0 0.1 0.0 0.0 0.1
## scale(not_pun_bonus) 0.0 0.1 -0.1 0.0 0.1
## sigma 1.0 0.0 1.0 1.0 1.1
##
## Fit Diagnostics:
## mean sd 10% 50% 90%
## mean_PPD 0.0 0.1 -0.1 0.0 0.1
##
## The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
##
## MCMC diagnostics
## mcse Rhat n_eff
## (Intercept) 0.0 1.0 4584
## scale(kappa) 0.0 1.0 3550
## scale(alpha_pos) 0.0 1.0 2929
## scale(alpha_neg) 0.0 1.0 4164
## scale(phi) 0.0 1.0 3188
## scale(rho) 0.0 1.0 3853
## scale(rl_off) 0.0 1.0 3976
## scale(epsilon) 0.0 1.0 4285
## scale(not_pun_bonus) 0.0 1.0 3743
## sigma 0.0 1.0 3779
## mean_PPD 0.0 1.0 4241
## log-posterior 0.1 1.0 1875
##
## For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).
packageVersion("rstanarm")
## [1] '2.32.1'
Get traces for each par
learn_traces <- data.frame(as.matrix(learn_preds_bayes))
learn_est <- ExtractEstimates(learn_traces)
learn_est$label <- "Learning"
learn_est$category <- "Proportion correct"
# Drop the coefs that aren't model pars
learn_est <- learn_est %>% filter(!coef %in% c("sigma", "X.Intercept."))
test_traces <- data.frame(as.matrix(test_preds_bayes))
test_est <- ExtractEstimates(test_traces)
test_est$label <- "Test"
test_est$category <- "Proportion correct"
test_est <- test_est %>% filter(!coef %in% c("sigma", "X.Intercept."))
assert(all(test_est$coef==learn_est$coef))
da_traces <- data.frame(as.matrix(da_preds_bayes))
da_est <- ExtractEstimates(da_traces)
da_est$label <- "Depr/anx"
da_est$category <- "Symptoms"
da_est <- da_est %>% filter(!coef %in% c("sigma", "X.Intercept."))
assert(all(da_est$coef==learn_est$coef))
ReturnPosteriorMeanAnd90CI(da_traces$scale.alpha_neg.)
##
## posterior mean = 0.005131044
## 90% CI =90% HDI: [-0.09, 0.11]
ReturnPosteriorMeanAnd90CI(da_traces$scale.alpha_pos.)
##
## posterior mean = 0.05431283
## 90% CI =90% HDI: [-0.06, 0.17]
ReturnPosteriorMeanAnd90CI(da_traces$scale.kappa.)
##
## posterior mean = -0.08965521
## 90% CI =90% HDI: [-0.21, 0.03]
ReturnPosteriorMeanAnd90CI(da_traces$scale.phi.)
##
## posterior mean = -0.07042382
## 90% CI =90% HDI: [-0.18, 0.05]
ReturnPosteriorMeanAnd90CI(da_traces$scale.rho.)
##
## posterior mean = 0.008827428
## 90% CI =90% HDI: [-0.10, 0.11]
ReturnPosteriorMeanAnd90CI(da_traces$scale.rl_off.)
##
## posterior mean = 0.005668033
## 90% CI =90% HDI: [-0.10, 0.12]
ReturnPosteriorMeanAnd90CI(da_traces$scale.epsilon.)
##
## posterior mean = 0.04158132
## 90% CI =90% HDI: [-0.06, 0.15]
ReturnPosteriorMeanAnd90CI(da_traces$scale.not_pun_bonus.)
##
## posterior mean = 0.06771418
## 90% CI =90% HDI: [-0.04, 0.17]
The three with close to some support
cat("\n Depression and bonus for non-punishing items\n")
##
## Depression and bonus for non-punishing items
ReturnPosteriorMeanAnd90CI(da_traces$scale.not_pun_bonus.)
##
## posterior mean = 0.06771418
## 90% CI =90% HDI: [-0.04, 0.17]
cat("\n Depression and WM decay\n")
##
## Depression and WM decay
ReturnPosteriorMeanAnd90CI(da_traces$scale.phi.)
##
## posterior mean = -0.07042382
## 90% CI =90% HDI: [-0.18, 0.05]
cat("\n Depression and WM capacity\n")
##
## Depression and WM capacity
ReturnPosteriorMeanAnd90CI(da_traces$scale.kappa.)
##
## posterior mean = -0.08965521
## 90% CI =90% HDI: [-0.21, 0.03]
cat("\n Depression and learning rate from negative PEs \n")
##
## Depression and learning rate from negative PEs
ReturnPosteriorMeanAnd90CI(da_traces$scale.alpha_neg.)
##
## posterior mean = 0.005131044
## 90% CI =90% HDI: [-0.09, 0.11]
Robustness check in univariate models the results aligned w hypotheses
summary(npb_da_preds_bayes <- rstanarm::stan_glm(
scale(depr_anx_sum) ~
scale(not_pun_bonus),
data=m35))
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 1.3e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.13 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.017 seconds (Warm-up)
## Chain 1: 0.034 seconds (Sampling)
## Chain 1: 0.051 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 6e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.016 seconds (Warm-up)
## Chain 2: 0.033 seconds (Sampling)
## Chain 2: 0.049 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 4e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.04 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.016 seconds (Warm-up)
## Chain 3: 0.033 seconds (Sampling)
## Chain 3: 0.049 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 6e-06 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.016 seconds (Warm-up)
## Chain 4: 0.032 seconds (Sampling)
## Chain 4: 0.048 seconds (Total)
## Chain 4:
##
## Model Info:
## function: stan_glm
## family: gaussian [identity]
## formula: scale(depr_anx_sum) ~ scale(not_pun_bonus)
## algorithm: sampling
## sample: 4000 (posterior sample size)
## priors: see help('prior_summary')
## observations: 275
## predictors: 2
##
## Estimates:
## mean sd 10% 50% 90%
## (Intercept) 0.0 0.1 -0.1 0.0 0.1
## scale(not_pun_bonus) 0.0 0.1 0.0 0.0 0.1
## sigma 1.0 0.0 1.0 1.0 1.1
##
## Fit Diagnostics:
## mean sd 10% 50% 90%
## mean_PPD 0.0 0.1 -0.1 0.0 0.1
##
## The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
##
## MCMC diagnostics
## mcse Rhat n_eff
## (Intercept) 0.0 1.0 3933
## scale(not_pun_bonus) 0.0 1.0 4136
## sigma 0.0 1.0 3796
## mean_PPD 0.0 1.0 3975
## log-posterior 0.0 1.0 1797
##
## For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).
Of note, results here and elsewhere very slightly discrepant when run within notebook vs. knitted results presumably due to rounding during knit/posterior sampling randomness but does not change results
cat("\n Univariate depression/anxiety bonus for non-punishing items\n")
##
## Univariate depression/anxiety bonus for non-punishing items
npb_traces <- data.frame(as.matrix(npb_da_preds_bayes))
ReturnPosteriorMeanAnd90CI(npb_traces$scale.not_pun_bonus.)
##
## posterior mean = 0.04539584
## 90% CI =90% HDI: [-0.06, 0.15]
rrs_traces <- data.frame(as.matrix(rrs_preds_bayes))
rrs_est <- ExtractEstimates(rrs_traces)
rrs_est$label <- "Rumination"
rrs_est$category <- "Symptoms"
rrs_est <- rrs_est %>% filter(!coef %in% c("sigma", "X.Intercept."))
assert(all(rrs_est$coef==test_est$coef))
ReturnPosteriorMeanAnd90CI(rrs_traces$scale.alpha_neg.)
##
## posterior mean = -0.002717084
## 90% CI =90% HDI: [-0.10, 0.11]
ReturnPosteriorMeanAnd90CI(rrs_traces$scale.alpha_pos.)
##
## posterior mean = 0.01449482
## 90% CI =90% HDI: [-0.09, 0.14]
ReturnPosteriorMeanAnd90CI(rrs_traces$scale.kappa.)
##
## posterior mean = -0.003328104
## 90% CI =90% HDI: [-0.13, 0.11]
ReturnPosteriorMeanAnd90CI(rrs_traces$scale.phi.)
##
## posterior mean = -0.02781697
## 90% CI =90% HDI: [-0.15, 0.09]
ReturnPosteriorMeanAnd90CI(rrs_traces$scale.rho.)
##
## posterior mean = 0.03524774
## 90% CI =90% HDI: [-0.07, 0.14]
ReturnPosteriorMeanAnd90CI(rrs_traces$scale.rl_off.)
##
## posterior mean = 0.04702406
## 90% CI =90% HDI: [-0.06, 0.16]
ReturnPosteriorMeanAnd90CI(rrs_traces$scale.epsilon.)
##
## posterior mean = 0.03527464
## 90% CI =90% HDI: [-0.07, 0.14]
ReturnPosteriorMeanAnd90CI(rrs_traces$scale.not_pun_bonus.)
##
## posterior mean = -0.001529133
## 90% CI =90% HDI: [-0.11, 0.10]
all_ests <- rbind(learn_est, test_est, da_est, rrs_est)
all_ests$label <- factor(all_ests$label, levels=c("Learning", "Test", "Depr/anx", "Rumination"))
all_ests_p <- ggplot(all_ests, #%>% filter(category=="Proportion correct"),
aes(y=coef, x=m, fill=label)) +
geom_vline(xintercept=0, size=1.3, color="gray35") +
geom_vline(xintercept=c(-.5, -.25, .25, .5), color="gray65") +
geom_errorbar(aes(xmin=lb_10, xmax=ub_90), width=.2, position=position_dodge(width=.4), size=1.5) +
geom_point(pch=21, size=5, alpha=.9, position=position_dodge(width=.4)) + ga + ap + lp +
ylab("Regression coefficient \n of model parameter") +
xlab("Mean and 90% credible interval") +
scale_y_discrete(labels=rev(c(
TeX('$\\RL^{off}'),
TeX('$\\rho'),
TeX('$\\phi'),
TeX('$\\n-p_{b}'),
TeX('$\\kappa'),
TeX('$\\epsilon'),
TeX('$\\alpha_{+}'),
TeX('$\\alpha_{-}')))) +
facet_wrap(~category) + ft +
scale_fill_manual(values=c("skyblue", "blue", "orange", "red")) +
scale_x_continuous(labels=function(x) sprintf("%.1f", x)) + theme(axis.text.y = element_text(size=32)) #+ tol
all_ests_p
Spot checks of plot
# hist(learn_traces$scale.kappa.)
# hist(test_traces$scale.kappa.)
# hist(da_traces$scale.kappa.)
# hist(learn_traces$scale.rl_off.)
# hist(test_traces$scale.rl_off.)
# hist(da_traces$scale.rl_off.)
# hist(test_traces$scale.alpha_pos.)
# hist(da_traces$scale.alpha_pos.)
# hist(learn_traces$scale.epsilon.)
# hist(test_traces$scale.epsilon.)
# ggsave("../paper/figs/pieces/fig8_modelpar_posteriors.png", all_ests_p ,
# width=14, height=12, dpi=500)
Robustness check in univariate models the results aligned w hypotheses
summary(phi_da_preds_bayes <- rstanarm::stan_glm(
scale(depr_anx_sum) ~
scale(phi),
data=m35))
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 1.2e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.12 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.018 seconds (Warm-up)
## Chain 1: 0.034 seconds (Sampling)
## Chain 1: 0.052 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 8e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.08 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.016 seconds (Warm-up)
## Chain 2: 0.034 seconds (Sampling)
## Chain 2: 0.05 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 7e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.07 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.015 seconds (Warm-up)
## Chain 3: 0.032 seconds (Sampling)
## Chain 3: 0.047 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 1e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.1 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.017 seconds (Warm-up)
## Chain 4: 0.033 seconds (Sampling)
## Chain 4: 0.05 seconds (Total)
## Chain 4:
##
## Model Info:
## function: stan_glm
## family: gaussian [identity]
## formula: scale(depr_anx_sum) ~ scale(phi)
## algorithm: sampling
## sample: 4000 (posterior sample size)
## priors: see help('prior_summary')
## observations: 275
## predictors: 2
##
## Estimates:
## mean sd 10% 50% 90%
## (Intercept) 0.0 0.1 -0.1 0.0 0.1
## scale(phi) 0.0 0.1 -0.1 0.0 0.0
## sigma 1.0 0.0 0.9 1.0 1.1
##
## Fit Diagnostics:
## mean sd 10% 50% 90%
## mean_PPD 0.0 0.1 -0.1 0.0 0.1
##
## The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
##
## MCMC diagnostics
## mcse Rhat n_eff
## (Intercept) 0.0 1.0 4126
## scale(phi) 0.0 1.0 4095
## sigma 0.0 1.0 3939
## mean_PPD 0.0 1.0 4023
## log-posterior 0.0 1.0 1756
##
## For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).
summary(kappa_da_preds_bayes <- rstanarm::stan_glm(
scale(depr_anx_sum) ~
scale(kappa),
data=m35))
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 1.1e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.11 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.017 seconds (Warm-up)
## Chain 1: 0.034 seconds (Sampling)
## Chain 1: 0.051 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 9e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.09 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.017 seconds (Warm-up)
## Chain 2: 0.032 seconds (Sampling)
## Chain 2: 0.049 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 8e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.08 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.017 seconds (Warm-up)
## Chain 3: 0.034 seconds (Sampling)
## Chain 3: 0.051 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 8e-06 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.08 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.017 seconds (Warm-up)
## Chain 4: 0.033 seconds (Sampling)
## Chain 4: 0.05 seconds (Total)
## Chain 4:
##
## Model Info:
## function: stan_glm
## family: gaussian [identity]
## formula: scale(depr_anx_sum) ~ scale(kappa)
## algorithm: sampling
## sample: 4000 (posterior sample size)
## priors: see help('prior_summary')
## observations: 275
## predictors: 2
##
## Estimates:
## mean sd 10% 50% 90%
## (Intercept) 0.0 0.1 -0.1 0.0 0.1
## scale(kappa) -0.1 0.1 -0.1 -0.1 0.0
## sigma 1.0 0.0 0.9 1.0 1.1
##
## Fit Diagnostics:
## mean sd 10% 50% 90%
## mean_PPD 0.0 0.1 -0.1 0.0 0.1
##
## The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
##
## MCMC diagnostics
## mcse Rhat n_eff
## (Intercept) 0.0 1.0 4011
## scale(kappa) 0.0 1.0 3801
## sigma 0.0 1.0 3926
## mean_PPD 0.0 1.0 3971
## log-posterior 0.0 1.0 1887
##
## For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).
summary(npb_da_preds_bayes <- rstanarm::stan_glm(
scale(depr_anx_sum) ~
scale(not_pun_bonus),
data=m35))
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 1e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.1 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.017 seconds (Warm-up)
## Chain 1: 0.032 seconds (Sampling)
## Chain 1: 0.049 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 7e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.07 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.016 seconds (Warm-up)
## Chain 2: 0.033 seconds (Sampling)
## Chain 2: 0.049 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 5e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.05 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.017 seconds (Warm-up)
## Chain 3: 0.033 seconds (Sampling)
## Chain 3: 0.05 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 3e-06 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.03 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.017 seconds (Warm-up)
## Chain 4: 0.033 seconds (Sampling)
## Chain 4: 0.05 seconds (Total)
## Chain 4:
##
## Model Info:
## function: stan_glm
## family: gaussian [identity]
## formula: scale(depr_anx_sum) ~ scale(not_pun_bonus)
## algorithm: sampling
## sample: 4000 (posterior sample size)
## priors: see help('prior_summary')
## observations: 275
## predictors: 2
##
## Estimates:
## mean sd 10% 50% 90%
## (Intercept) 0.0 0.1 -0.1 0.0 0.1
## scale(not_pun_bonus) 0.0 0.1 0.0 0.0 0.1
## sigma 1.0 0.0 0.9 1.0 1.1
##
## Fit Diagnostics:
## mean sd 10% 50% 90%
## mean_PPD 0.0 0.1 -0.1 0.0 0.1
##
## The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
##
## MCMC diagnostics
## mcse Rhat n_eff
## (Intercept) 0.0 1.0 3878
## scale(not_pun_bonus) 0.0 1.0 4073
## sigma 0.0 1.0 4175
## mean_PPD 0.0 1.0 4030
## log-posterior 0.0 1.0 1984
##
## For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).
Sanity check against predictions of frequentist models and check their VIF
summary(da_preds <- lm(scale(depr_anx_sum) ~
scale(kappa) +
scale(alpha_pos) +
scale(alpha_neg) +
scale(phi) +
scale(rho) +
scale(rl_off) +
scale(epsilon) +
scale(not_pun_bonus),
data=m35))
##
## Call:
## lm(formula = scale(depr_anx_sum) ~ scale(kappa) + scale(alpha_pos) +
## scale(alpha_neg) + scale(phi) + scale(rho) + scale(rl_off) +
## scale(epsilon) + scale(not_pun_bonus), data = m35)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.2736 -0.8104 -0.2612 0.7150 3.0819
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.124e-16 6.069e-02 0.000 1.000
## scale(kappa) -9.197e-02 7.175e-02 -1.282 0.201
## scale(alpha_pos) 5.706e-02 7.000e-02 0.815 0.416
## scale(alpha_neg) 4.365e-03 6.299e-02 0.069 0.945
## scale(phi) -7.165e-02 7.138e-02 -1.004 0.316
## scale(rho) 6.362e-03 6.506e-02 0.098 0.922
## scale(rl_off) 3.362e-03 6.608e-02 0.051 0.959
## scale(epsilon) 4.345e-02 6.302e-02 0.689 0.491
## scale(not_pun_bonus) 6.732e-02 6.285e-02 1.071 0.285
##
## Residual standard error: 1.006 on 266 degrees of freedom
## Multiple R-squared: 0.01666, Adjusted R-squared: -0.01292
## F-statistic: 0.5633 on 8 and 266 DF, p-value: 0.8076
summary(rrs_preds <- lm(scale(RRS_sum) ~
scale(kappa) +
scale(alpha_pos) +
scale(alpha_neg) +
scale(phi) +
scale(rho) +
scale(rl_off) +
scale(epsilon) +
scale(not_pun_bonus),
data=m35))
##
## Call:
## lm(formula = scale(RRS_sum) ~ scale(kappa) + scale(alpha_pos) +
## scale(alpha_neg) + scale(phi) + scale(rho) + scale(rl_off) +
## scale(epsilon) + scale(not_pun_bonus), data = m35)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.54304 -0.80805 -0.03143 0.61693 2.77658
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -7.697e-17 6.099e-02 0.000 1.000
## scale(kappa) -2.798e-03 7.211e-02 -0.039 0.969
## scale(alpha_pos) 1.444e-02 7.035e-02 0.205 0.838
## scale(alpha_neg) -2.077e-03 6.330e-02 -0.033 0.974
## scale(phi) -2.554e-02 7.173e-02 -0.356 0.722
## scale(rho) 3.520e-02 6.538e-02 0.538 0.591
## scale(rl_off) 4.736e-02 6.641e-02 0.713 0.476
## scale(epsilon) 3.551e-02 6.333e-02 0.561 0.575
## scale(not_pun_bonus) -3.289e-03 6.316e-02 -0.052 0.959
##
## Residual standard error: 1.011 on 266 degrees of freedom
## Multiple R-squared: 0.006906, Adjusted R-squared: -0.02296
## F-statistic: 0.2312 on 8 and 266 DF, p-value: 0.9849
summary(learning_preds <- lm(scale(learn_m) ~
scale(kappa) +
scale(alpha_pos) +
scale(alpha_neg) +
scale(phi) +
scale(rho) +
scale(rl_off) +
scale(epsilon) +
scale(not_pun_bonus),
data=m35))
##
## Call:
## lm(formula = scale(learn_m) ~ scale(kappa) + scale(alpha_pos) +
## scale(alpha_neg) + scale(phi) + scale(rho) + scale(rl_off) +
## scale(epsilon) + scale(not_pun_bonus), data = m35)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.94486 -0.24830 0.04602 0.28117 1.60594
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.684e-15 3.183e-02 0.000 1.000000
## scale(kappa) 2.549e-01 3.764e-02 6.772 8.11e-11 ***
## scale(alpha_pos) 4.301e-01 3.672e-02 11.715 < 2e-16 ***
## scale(alpha_neg) 1.557e-01 3.304e-02 4.713 3.93e-06 ***
## scale(phi) -2.687e-01 3.744e-02 -7.176 7.13e-12 ***
## scale(rho) 1.363e-01 3.412e-02 3.995 8.38e-05 ***
## scale(rl_off) -2.667e-01 3.466e-02 -7.694 2.78e-13 ***
## scale(epsilon) -3.080e-01 3.306e-02 -9.319 < 2e-16 ***
## scale(not_pun_bonus) -1.135e-01 3.297e-02 -3.444 0.000666 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5279 on 266 degrees of freedom
## Multiple R-squared: 0.7295, Adjusted R-squared: 0.7213
## F-statistic: 89.65 on 8 and 266 DF, p-value: < 2.2e-16
summary(testing_preds <- lm(scale(test_m) ~
scale(kappa) +
scale(alpha_pos) +
scale(alpha_neg) +
scale(phi) +
scale(rho) +
scale(rl_off) +
scale(epsilon) +
scale(not_pun_bonus),
data=m35))
##
## Call:
## lm(formula = scale(test_m) ~ scale(kappa) + scale(alpha_pos) +
## scale(alpha_neg) + scale(phi) + scale(rho) + scale(rl_off) +
## scale(epsilon) + scale(not_pun_bonus), data = m35)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.0490 -0.2169 0.0521 0.3238 1.3265
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.910e-15 3.284e-02 0.000 1.0000
## scale(kappa) 5.739e-02 3.883e-02 1.478 0.1406
## scale(alpha_pos) 6.386e-01 3.788e-02 16.858 < 2e-16 ***
## scale(alpha_neg) 1.718e-01 3.408e-02 5.040 8.63e-07 ***
## scale(phi) -1.834e-01 3.862e-02 -4.747 3.37e-06 ***
## scale(rho) 7.813e-02 3.520e-02 2.219 0.0273 *
## scale(rl_off) -4.875e-01 3.576e-02 -13.633 < 2e-16 ***
## scale(epsilon) -1.670e-01 3.410e-02 -4.896 1.70e-06 ***
## scale(not_pun_bonus) -3.211e-02 3.401e-02 -0.944 0.3460
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5446 on 266 degrees of freedom
## Multiple R-squared: 0.7121, Adjusted R-squared: 0.7034
## F-statistic: 82.22 on 8 and 266 DF, p-value: < 2.2e-16
car::vif(learning_preds)
## scale(kappa) scale(alpha_pos) scale(alpha_neg)
## 1.392691 1.325469 1.073142
## scale(phi) scale(rho) scale(rl_off)
## 1.378196 1.144827 1.181171
## scale(epsilon) scale(not_pun_bonus)
## 1.074371 1.068519
car::vif(testing_preds)
## scale(kappa) scale(alpha_pos) scale(alpha_neg)
## 1.392691 1.325469 1.073142
## scale(phi) scale(rho) scale(rl_off)
## 1.378196 1.144827 1.181171
## scale(epsilon) scale(not_pun_bonus)
## 1.074371 1.068519
car::vif(da_preds)
## scale(kappa) scale(alpha_pos) scale(alpha_neg)
## 1.392691 1.325469 1.073142
## scale(phi) scale(rho) scale(rl_off)
## 1.378196 1.144827 1.181171
## scale(epsilon) scale(not_pun_bonus)
## 1.074371 1.068519
car::vif(rrs_preds)
## scale(kappa) scale(alpha_pos) scale(alpha_neg)
## 1.392691 1.325469 1.073142
## scale(phi) scale(rho) scale(rl_off)
## 1.378196 1.144827 1.181171
## scale(epsilon) scale(not_pun_bonus)
## 1.074371 1.068519
m41_v1 <- read.csv("../model_res/opt/best/BEST__RunHWMPRew34387.csv")
m41_v2 <- read.csv("../model_res/opt/best/BEST__RunHWMPRew47450.csv")
m41 <- rbind(m41_v1, m41_v2) %>% group_by(ID) %>% slice(which.min(nll))
#write.csv(m41, "../model_res/opt/best/BEST__m41_RunHWMPRew.csv")
hist(m41$learning_bias, breaks=100)
hist(m41$alpha_ck, breaks=100)
hist(m35$alpha_pos, breaks=100)
Highly correlated
ComparePars(m35$alpha_pos, m41$alpha_ck)
As expected not very correlated
ComparePars(m35$alpha_neg, m41$alpha_ck)
hist(m41$ck_off, breaks=100)
median(m41$ck_off)
## [1] 0.07782407
median(m35$rl_off)
## [1] 0.6136962
ComparePars(m35$rl_off, m41$ck_off)
Same number of pars so don’t need penalty to compare
ComparePars(m35$nll, m41$nll, "", "m35", "m41")
sum(m35$nll-m41$nll)
## [1] -500.0767
length(which(m35$nll < m41$nll))/length(m41$nll)
## [1] 0.6581818
0 rather than 1/3 inits
m42_v1 <- read.csv("../model_res/opt/best/BEST__RunHWMPRew0Inits18746.csv")
m42_v2 <- read.csv("../model_res/opt/best/BEST__RunHWMPRew0Inits77616.csv")
m42 <- rbind(m42_v1, m42_v2) %>% group_by(ID) %>% slice(which.min(nll))
Very similar but slightly worse than m41 and thus still substantially worse than m35
ComparePars(m35$nll, m42$nll, "", "m35", "m42")
sum(m35$nll-m42$nll)
## [1] -510.909
length(which(m35$nll < m42$nll))/length(m42$nll)
## [1] 0.6618182
m41_s <- read.csv("../model_res/sims/SIM_RunHWMPRew25874.csv")
m41_s_learn <- m41_s %>% filter(type=="learning")
m41_s_test <- m41_s %>% filter(type=="test")
pcor_ss_sim_m41 <- data.frame(m41_s_learn %>% group_by(set_size, stim_iter) %>%
summarize(m=mean(corrects), n()))
## `summarise()` has grouped output by 'set_size'. You can override using the
## `.groups` argument.
pcor_ss_sim_m41_iter <- data.frame(m41_s_learn %>% filter(iter %in% c(1:30)) %>% group_by(stim_iter, set_size, iter) %>%
summarize(m=mean(corrects), n()))
## `summarise()` has grouped output by 'stim_iter', 'set_size'. You can override
## using the `.groups` argument.
sim_m41_p1 <-
ggplot(pcor_ss_sim_m41, aes(x=as.factor(stim_iter), y=m, group=as.factor(set_size), color=as.factor(set_size))) +
geom_line() +
geom_hline(yintercept=.33, size=1.5, color="gray57") + # chance line
geom_hline(yintercept=c(.5, .6, .7, .8, .9, 1), linetype="dotted") +
geom_vline(xintercept=c(2, 5, 8, 10), linetype="dotted") +
geom_jitter(data=pcor_ss_sim_m41_iter, aes(fill=as.factor(set_size)), color="black", height=0, width=.2, alpha=1, size=2, pch=21) +
geom_point(aes(fill=as.factor(set_size)), color="black", size=6, pch=21, alpha=.7) +
annotate("rect", xmin=6, xmax=10.5, ymin=.3, ymax=1.1, alpha=0.2, fill="gray57") +
ga + ap + lp + xlab("Stimulus iteration") + ylab("Proportion correct") +
tp + ggtitle("Simulated")
#ggsave("../paper/figs/pieces/supp6_sim_perf.png", sim_m41_p1, height = 3.5, width=6, dpi=300)
emp_sim_perf <- emp_p1 + sim_m41_p1
emp_sim_perf
np_si_sims_summs <- m41_s %>%
filter(type=="learning" & corrects==0) %>% group_by(stim_iter) %>% #filter(stim_iter %in% c(2:10)) %>%
summarize(mw=mean(worsts), mn=mean(neutrals), n=n()) %>% mutate(neutral_pref=(mn-mw))
np_si_sims_summs_var <- m41_s %>%
filter(type=="learning" & corrects==0) %>% group_by(stim_iter, iter) %>% #filter(stim_iter %in% c(2:10)) %>%
summarize(mw=mean(worsts), mn=mean(neutrals), n=n()) %>% mutate(neutral_pref=(mn-mw))
## `summarise()` has grouped output by 'stim_iter'. You can override using the
## `.groups` argument.
m41_np_si_plot <- ggplot(np_si_sims_summs,
aes(x=stim_iter, y=neutral_pref, fill=as.numeric(stim_iter))) +
geom_hline(yintercept=.0, size=1.5, color="gray57") + # chance line
geom_bar(stat="identity", color="black") +
geom_jitter(data=np_si_sims_summs_var, aes(x=stim_iter, y=neutral_pref), height=0,
width=.2, pch=21) +
annotate("rect", xmin=5.5, xmax=10.5, ymin=0, ymax=.22, alpha=0.2, pch=21) +
annotate("text", x=3, y=.20, label="", size=8) +
annotate("text", x=7.5, y=.20, label="", size=8) +
ga + ap + tol + xlab("Simulus iteration") + ylab("") + tp +
ylim(-.05, .23) +
ylab("Neutral preference \n during learning") +
ggtitle("Simulated") + scale_color_gradient2() +
scale_x_continuous(breaks=seq(1, 10, 1)) + ylim(-.1, .23)
## Warning in annotate("rect", xmin = 5.5, xmax = 10.5, ymin = 0, ymax = 0.22, :
## Ignoring unknown parameters: `shape`
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
m41_np_si_plot
#ggsave("../paper/figs/pieces/supp6_neut-pref_emp-and-sims.png", m41_np_si_plot, height = 4, width=11, dpi=300)
Test and SI6 data
pcor_p1_test_sim_m41 <-
data.frame(m41_s_test %>% filter(phase==1) %>% group_by(set_size) %>% summarize(m=mean(corrects)))
pcor_p1_test_sim_m41_iters <-
data.frame(m41_s_test %>% filter(phase==1) %>% group_by(set_size, iter) %>% summarize(m=mean(corrects)))
## `summarise()` has grouped output by 'set_size'. You can override using the
## `.groups` argument.
pcor_si6_sim_m41 <-
data.frame(m41_s_learn %>% filter(stim_iter==6) %>% group_by(set_size) %>% summarize(m=mean(corrects)))
pcor_si6_sim_m41_iters <-
data.frame(m41_s_learn %>% filter(stim_iter==6) %>% group_by(set_size, iter) %>% summarize(m=mean(corrects)))
## `summarise()` has grouped output by 'set_size'. You can override using the
## `.groups` argument.
pcor_p2_test_sim_m41 <-
data.frame(m41_s_test %>% filter(phase==2) %>% group_by(set_size) %>% summarize(m=mean(corrects)))
pcor_p2_test_sim_m41_iters <-
data.frame(m41_s_test %>% filter(phase==2) %>% group_by(set_size, iter) %>% summarize(m=mean(corrects)))
## `summarise()` has grouped output by 'set_size'. You can override using the
## `.groups` argument.
sim_p1_test_m41 <-
ggplot(pcor_p1_test_sim_m41, aes(x=as.factor(set_size), y=m, fill=as.factor(set_size))) +
geom_hline(yintercept = seq(.1, 1, .1), alpha=.3) +
geom_bar(stat="identity", color="black") +
geom_jitter(data=pcor_p1_test_sim_m41_iters, size=2, alpha=1, width=.08, height=0, pch=21,
aes(x=as.factor(set_size), y=m, fill=as.factor(set_size))) +
ga + ap + xlab("Set size") + ylab("Proportion correct") + tol + ylim(0, 1) +
tp + ggtitle("Test phase 1")
sim_si6_m41 <- ggplot(pcor_si6_sim_m41, aes(x=as.factor(set_size), y=m, fill=as.factor(set_size))) +
geom_hline(yintercept = seq(.1, 1, .1), alpha=.3) +
geom_bar(stat="identity", color="black") +
geom_jitter(data=pcor_si6_sim_m41_iters, size=2, alpha=1, width=.08, height=0, pch=21,
aes(x=as.factor(set_size), y=m, fill=as.factor(set_size))) +
ga + ap + xlab("Set size") + ylab("") + tol + ylim(0, 1) +
tp + ggtitle("Stimulus iteration 6")
sim_p2_test_m41 <- ggplot(pcor_p2_test_sim_m41, aes(x=as.factor(set_size), y=m, fill=as.factor(set_size))) +
geom_hline(yintercept = seq(.1, 1, .1), alpha=.3) +
geom_bar(stat="identity", color="black") +
geom_jitter(data=pcor_p2_test_sim_m41_iters, size=2, alpha=1, width=.08, height=0, pch=21,
aes(x=as.factor(set_size), y=m, fill=as.factor(set_size))) +
ga + ap + xlab("Set size") + ylab("") + tol + ylim(0, 1) + tp +
tp + ggtitle("Test phase 2")
sims_u_plot <-
sim_p1_test_m41 + sim_si6_m41 + sim_p2_test_m41 + plot_annotation(title="Simulated", theme = theme(plot.title = element_text(size = 25, hjust=.5)))#,
emps_u_plot
sims_u_plot
#ggsave("../paper/figs/pieces/supp6_simU.png", sims_u_plot, height = 4, width=11, dpi=300)
hist(m35$rl_off, breaks=50)
ck_off_p <-
ggplot(m41, aes(x=ck_off)) +
geom_vline(xintercept=median(m41$ck_off), size=4) +
geom_histogram(fill="white", color="black") + ga + ap + ylab("") +
ggtitle(TeX('$\\CK^{off}')) + tp + xlab("")
ck_off_p
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
median(m41$ck_off)
## [1] 0.07782407
median(m35$rl_off)
## [1] 0.6136962
Note: the learning bias par was implemented as a scalar on WM and RL such that lower values correspond to more shrinking of parameter, so plotting 1-bias so it has more intuitive interpretation as higher value = more insensitivity of WM and CK after negative PEs
lb_off_p <-
ggplot(m41, aes(x=1-learning_bias)) +
geom_vline(xintercept=median(1-m41$learning_bias), size=4) +
geom_histogram(fill="white", color="black") + ga + ap + ylab("") +
ggtitle(TeX('bias')) + tp + xlab("")
lb_off_p
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Test neutral pref
m41_test_sim_error_df <- m41_s_test %>% filter(corrects==0)
np_test_summs_m41_sim_var <- m41_test_sim_error_df %>% group_by(phase, iter) %>%
summarize(mw=mean(worsts), mn=mean(neutrals), n=n()) %>% mutate(neutral_pref=(mn-mw))
## `summarise()` has grouped output by 'phase'. You can override using the
## `.groups` argument.
np_test_summs_m41_sim <- m41_test_sim_error_df %>% group_by(phase) %>%
summarize(mw=mean(worsts), mn=mean(neutrals), n=n()) %>% mutate(neutral_pref=(mn-mw))
sim_m41_np_si_test_plot <-
ggplot(np_test_summs_m41_sim , aes(x=as.factor(phase), y=neutral_pref, fill=as.factor(phase))) +
geom_hline(yintercept=.0, size=1.5, color="gray57") + # chance line
geom_jitter(data=np_test_summs_m41_sim_var, height=0, width=.2, size=3, alpha=.8) +
geom_bar(stat="identity", color="black", alpha=.8) +
ga + ap + tol + xlab("Phase") + ylab("") + tp +
ylab("") +
ggtitle("Neutral Preference at Test") + scale_fill_manual(values=c("gray81", "gray40")) + ylim(-.05, .23)
sim_m41_np_si_test_plot
#ggsave("../paper/figs/pieces/supp6_emp_np_si.png", sim_m41_np_si_test_plot, height = 4, width=5, dpi=300)